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§0.  Introduction. 

The  objective  of  this  research  is  to  provide  new  analytic  approaches  and  new  mathe¬ 
matical  tools  for  combat  modeling  and  simulation.  In  this  paper  we  focus  on  three  methods 
for  assigning  values  (scores,  importances)  to  various  weapons  in  a  simulated  engagement. 
The  potential/antipotential  or  eigenvalue  method  and  the  method  used  by  the  U.S.  Army 
Concepts  Analysis  Agency  as  part  of  its  ATCAL  (Attrition  Calibration)  procedure,  are 
part  of  current  practice.  The  fire  allocation  methods  are,  relatively  speaking,  a  new  idea 
and  lead  to  some  very  interesting  systems  of  non-linear  equations  that  can  be  viewed  as 
non-linear  generalizations  of  the  simultaneous  eigenvalue  problem  that  occurs  in  the  (lin¬ 
ear)  eigenvalue  method.  This  similarity  with  the  eigenvalue  method  is  striking  and  leads 
one  to  conclude  that  fire  allocation  methods  are  the  “right”  non-linear  generalization  of 
the  eigenvalue  method. 

We  have  been  able  to  derive  a  number  of  useful  results.  For  example,  as  part  of  our 
investigation  of  the  ATCAL  weapon  importances,  we  were  able  to  establish  the  existence 
and  uniqueness  of  the  solution  to  the  importance  equations.  The  proof,  which  is  amazingly 
similar  to  the  proof  of  the  Perron- Frobenius  Theorem  that  gives  existence  and  uniqueness 
for  the  eigenvalue  method,  provides  a  new  iterative  scheme  for  solving  the  ATCAL  impor¬ 
tance  equations.  Again,  the  similarity  with  aspects  of  the  (linear)  eigenvalue  method  is 
interesting,  and  one  is  led  to  conclude  that  the  ATCAL  importances  will  inherit  many  of 
the  same  problems  that  occur  with  the  eigenvalue  method. 

We  also  discuss  various  methods  for  dealing  with  the  systems  of  non-linear  equations 
that  occur  in  the  fire  allocation  methods.  These  problems  take  us  to  the  frontiers  of  modern 
mathematical  research  into  computational  methods  for  finding  resultants  and  Grobner 
bases. 

Finally,  we  conclude  with  some  numerical  experiments  that  tested  the  sensitivity  of 
the  three  scoring  methods  to  changes  in  various  inputs  and  that  made  some  comparisons 
between  them.  Because  large  examples  would  be  computationally  expensive,  these  exper¬ 
iments  were  run  only  for  small  examples  having  at  most  three  weapons  per  side. 

Clearly  warfare  involves  many  dynamically  interacting  “variables”  and  a  large  mea¬ 
sure  of  unpredictable  human  behaviour.  As  a  result,  it  is  not  all  that  suited  to  scien- 
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tific/mathematical  analysis.  Nevertheless,  we  feel  that  it  is  worthwhile  to  examine  a  variety 
of  mathematical  paradigms  and  techniques,  so  as  to  be  able  to  build  the  richest  possible 
models  that  can  accommodate  the  widest  range  of  phenomena.  Moreover,  the  mathemat¬ 
ical  imderstanding  of  a  model’s  behavior  is  critical  to  effectively  moving  results  between 
models  in  a  hierarchy.  For  those  reasons,  basic  research  into  foundational  issues  in  combat 
modeling  and  simulation  is  important. 

The  author  would  like  to  thank  the  U.S.  Army  Concepts  Analysis  Agency,  particularly 
the  Chief  of  its  Advanced  Research  Projects  Office  (ARPO),  Mr.  Gerald  E.  Cooper,  for 
their  support  of  this  research  and  their  hospitality  during  visits  to  the  agency. 
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§1.  An  Analysis  of  the  Potential/ Anti-Potential  or  Eigenvalue  Method  for 
Assigning  Values  (Scores)  to  Weapon  Types. 

a.  The  Method. 

Numerous  methods  have  been  proposed  for  assigning  “values”  to  weapons  in  combat 
models  (Anderson  [56]).  Typically  the  resources  of  the  two  opposing  sides  are  grouped 
into  a  number  of  generic  types  (e.g.  MlAl,  MLRS,  BMP,  TRUCK,  etc.),  with  the  object 
of  the  method  being  to  compute  a  value  (or  score)  for  each  weapon  of  a  particular  type. 
These  values  are  then  used  to  make  comparisons  across  types  and  to  make  overall  force 
comparisons. 

Our  purpose  in  considering  here  the  Potential/ Anti-Potential  method  and  some  of  its 
variants  is  three- fold.  First,  we  will  use  it  to  illustrate  how  some  rather  sophisticated  math¬ 
ematical  ideas  and  important  computational  issues  can  arise  in  a  relatively  simple  context. 
For  example,  Anderson  [56]  points  out  that  some  methods  of  assigning  weapons  scores  are 
in  everyday  use,  eventhough  it  is  not  known  whether  solutions  of  the  equations  used  in  that 
method  always  exist,  and  if  they  exist,  whether  or  not  they  are  unique  (e.g.  ATCAL,  see 
Section  3  below).  Second,  it  is  always  important  to  know  the  limitations  of  your  mathemat¬ 
ical  methodology.  Fortunately,  with  the  Potential/ Anti-Potential  method  we  can  describe 
those  limitations  precisely.  In  addition,  whenever  one  selects  a  particular  mathematical 
paradigm  to  use  in  building  a  model,  one  frequently  can  derive  additional  properties  of 
the  model  as  consequences  of  the  initial  assumptions.  These  derived  features  can  serve  as 
a  kind  of  check  on  the  reasonableness  of  the  choice  of  mathematical  technique/paradigm. 
Such  “derived  consequences”  can  be  illustrated  with  the  Potential/ Anti-potential  method. 
Our  final  purpose  in  considering  this  method  is  that  many  of  its  features  carry  over  to  the 
other  two  methods  we  will  be  discussing  —  fire  allocation  and  ATCAL  importances.  Thus 
this  method  will  serve  to  illustrate  some  of  our  later  ideas.  For  all  of  these  reasons,  we 
have  chosen  to  discuss  the  method  in  considerable  detail. 

We  suppose  that  the  two  opposing  forces  (Blue  and  Red)  are  made  up  of  a  mix  of 
weapons  systems  broken  down  into  m  types  for  the  Blue  force  and  n  types  for  the  Red 
force.  Let  Bi  denote  the  number  of  weapons  of  type  i  in  the  Blue  force  (i  =  1, . . . ,  m)  and 
let  Rj  denote  the  number  of  weapons  of  type  j  in  the  Red  force  {j  =  1, . . . ,  n).  We  denote 
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the  sought  after  values  by  VB,  and  VRj.  Thus  VB,  (V Rj)  represents  the  value  per  type 
i  (;■)  weapon  in  the  Blue  (Red)  force.  We  also  let  Kij  denote  the  rate  at  which  one  Blue 
weapon  of  type  i  kills  Red  weapons  of  type  j  and  let  Lji  denote  the  rate  at  which  one  Red 
weapon  of  type  j  kills  Blue  weapons  of  type  i.  Kij  is  therefore  measured  in  units  of  Red 
weapons  of  type  j  per  Blue  weapon  of  type  i  per  unit  of  time. 

In  a  very  general  setting,  the  Kij  and  Lji  would  be  functions  of  the  Bi  s  and  Rj  s  so 
that  the  values  VBi  and  VRj  would  change  as  the  composition  of  the  Red  and  Blue  forces 
changed.  To  the  extent  that  the  Kij  and  Lji  are  non-constant,  the  total  value  for  each 
force 

B  =  Y^  VBi  xBi  and  R  =  J2  ^ 
i=l  i=l 

would  be  a  non-linear  function  of  the  force  components  Bi  and  Rj  and  of  time.  However 
the  eigenvalue  method  can  be  viewed  as  instantaneous  so  that  Kij  and  Lji  may  be  taken 
as  constant. 

The  basic  assumption  of  the  potential/anti-potential  or  eigenvalue  method  is  that  the 
Vcilue  per  weapon  of  a  particular  type  is  proportional  to  the  rate  at  which  a  weapon  of  that 
type  destroys  the  value  of  the  enemy’s  weapons.  We  can  assume  a  single  positive  constant 
of  proportionality  C  or  separate  positive  constants  CB  and  CR  for  the  Blue  and  Red  forces 
respectively.  We  also  impose  the  obvious  constraint  that  all  the  values  be  non-negative. 
This  leads  to  the  following  system  of  m  -|-  n  equations: 

n 

CB  xVBi=J2  ^<ijVRj  z  =  1, . . . ,  m 

i=i 

m 

CR  X  VRj  =  Y,  LjiVRi  i  =  1,  • .  • ,  n 

i=l 


or  in  matrix  terms: 


(1)  CB  xVB  =  K  X  VR 

(2)  CRxVR  =  LxVB 

subject  to  the  constraints  VBi  >  0  for  i  =  1, ...  ,Tn,  VRj  >  0  for  j  =  1, . . . , n,  CR  >  0, 
and  CB  >  0.  Note  that  all  the  Kij  and  Lji  are  always  non-negative  and  that  CB  and 
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CR,  in  addition  to  VB  and  VR,  are  unknowns.  Substituting  (2)  into  (1)  and  vice  versa 
yields  the  simultaneous  eigenvalue  problem 

(3)  \VB  =  Mkl  VB 

(4)  \VR^  Mlk  VR 

where  A  =  CRy.  CB  is  a  scalar,  Mkl  is  the  m  x  m  matrix  product  of  K  with  L,  and  Mlk 
is  the  n  X  n  matrix  product  of  L  with  K.  Note  that  both  Mkl  and  Mlk  are  non-negative, 
that  is,  have  non-negative  entries.  For  future  reference,  we  also  note  that  this  system 
(3),  (4)  is  remarkably  similar  to  a  system  of  non-linear  equations  that  occurs  in  the  fire 
allocation  method  of  Section  2. 

Because  Mkl  and  Mlk  are  obtained  as  the  product  of  K  and  L  in  different  orders, 
it  can  be  shown  that  they  have  the  same  non- zero  eigenvalues  (counted  with  multiplicity). 
Moreover,  the  Perron-Frobenius  Theorem  (see  below)  implies  that  if  Mkl  and  Mlk  are 
irreducible,  then  our  system  (3)  and  (4)  will  have  an  essentially  unique  solution.  Specif¬ 
ically,  Mkl  and  Mlk  will  have  a  common  positive  real  eigenvalue  A  whose  associated 
eigenvectors  VB  and  VR  will  have  all  positive  entries.  This  eigenvalue  A  =  Amax  will  be 
the  largest  real  eigenvalue  of  both  Mkl  and  Mlk-,  and  the  corresponding  eigenvectors 
VB  and  VR  will  be  unique  up  to  a  positive  scale  factor.  Simply  put,  there  will  be  (up  to 
scale)  only  one  eigenvector  for  each  of  Mkl  and  Mlk  with  non-negative  entries,  and  both 
will  correspond  to  the  largest  eigenvalue.  If  either  matrix  is  reducible,  then  zero  values  are 
possible,  and  the  solution  may  no  longer  be  unique. 

It  is  important  to  note  that  any  choice  of  positive  eigenvectors  V B  and  V R  will 
provide  a  solution  to  our  original  system  (1)  and  (2).  To  see  this,  suppose  we  have  selected 
two  such  positive  eigenvectors  V B  and  V R,  then 

Mkl  ^(Kx  VR)  =  {K  x  L)  x  (K  x  VR)  =  K  x  MlkVR  =  X{K  x  VR) 

so  that  K  X  VR  is  also  a  non-zero  non-negative  eigenvector  of  Mkl-  By  the  uniqueness 
result  mentioned  above,  K  xVR  must  therefore  be  a  positive  multiple,  call  it  CB,  of  VB, 
so  that 

CB  xVB  =  K  xVR 
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as  desired.  A  similar  calculation  gives 


CRxVR  =  L  X  VB. 

If  instead,  we  select  VR  =  r  VR,  r  >  0  and  VB  =  b  VB,  b  >  0  as  out  eigenvectors,  then 

^  xVB  =  K  xVR 

and 

CRxVR  =  LxVB 

where  ^  f  and  ^  =  CiZ  •  so  that  ^  CB=^CR  CB  =  A^ax  remains 

invariant.  Notice  that  if  6  =  r  then  CR  and  CB  do  not  change.  We  thus  have  two  ways 
to  view  our  two  degrees  of  freedom:  we  can  freely  scale  the  eigenvectors  V B  and  V R 
independently,  or  we  can  scale  CB  freely  (and  hence  CA  reciprocally)  and  then  scale  VB 
freely  (in  which  event  VR  is  uniquely  determined). 

Thus  all  we  really  have  at  this  point  is  relative  weapon  values  for  Blue  and  Red 
^  and  with  no  way  to  compare  Blue  weapons  to  Red  weapons.  One  common 

normalization  which  overcomes  this,  is  to  set  CB  —  CR  —  y/Kalx  (which  in  effect  equates 
units  of  Red  value  with  units  of  Blue  value).  One  can  then  do  relative  comparisons  across 
Red  and  Blue  weapons  types.  Finally,  to  obtain  fixed  values,  one  can  normalize  one  weapon 
type  to  have  value  =  1,  which  in  effect  amounts  to  a  change  of  scale  in  the  units  of  “value”. 
Whether  we  chose  to  normalize  or  not,  it  is  easy  to  see  that  the  quantity 

^/CRxR 
VCB  X  B' 

n  n 

where  R  =  V Rj  x  Rj  and  B  =  V Bi  x  Bi  are  the  total  Red  and  Blue  values 

j=i  t=i 

respectively,  is  constant  independent  of  our  two  scaling  degrees  of  freedom.  We  shall  see 
below  that  this  quantity  has  a  natural  interpretation  as  the  effective  force  ratio.  It  is  a 
combination  of  the  quantitative  ratio  §  and  the  qualitative  ratio  Note  that  doubling 
quantitative  strength  is  equivalent  to  a  four-fold  increase  in  qualitative  strength. 

We  remark  that  other  normalization  schemes  have  been  criticized  (Anderson  [56])  on 
the  grounds  that  dividing  a  weapon  type  into  two  “very  similar”  types  should  not  upset 
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the  force  ratio.  Those  arguments  appear  to  be  based  on  interpreting  ^  as  the  force  ratio, 
which  is  not  the  appropriate  measure.  Finally,  note  that  redefining  value  to  be  \/CR  VRj 
and  y/C BVBi  (i.e.  r  =  y/CR  and  b  =  y/CB)  is  equivalent  to  setting  CR  =  CB  =  VAj^ax- 
In  general  CB{CR)  involves  units  of  (time)-^  and  can  be  regarded  as  the  rate  at 
which  one  unit  of  Blue  (Red)  value  kills  units  of  Red  (Blue)  value.  This  is  because 

T.KiiVRj 

CB  =  ^  -  for  every  i  = 

V  IJ  i 


so  that 


ZBiZIUjxVRj 
CB  =  —  ^ 


Y.B^xVB, 

i 

The  numerator  is  the  units  of  Red  value  killed  per  unit  time,  and  the  denominator  is  total 
Blue  value.  Thus  CB  is  the  number  of  units  of  Red  value  killed  per  unit  of  Blue  value 
per  unit  time.  Hence,  ii  CR  and  CB  are  different,  the  units  of  Red  and  Blue  value  are 
qualitatively  different. 

The  invariant  quantity  y/CB  x  CR  =  VAmax  =  I  can  be  regarded  as  a  measure  of 
the  intensity  of  combat.  If  we  scale  all  kill  rates  Kij  and  Lij  by  c,  then  K  x  L  and  L  x  K 
scale  by  c^,  and  the  common  maximum  eigenvalue  changes  by  a  factor  of  c^.  This  in  turn 
means  that  I  scales  by  c,  which  is  appropriate  for  a  measure  of  intensity.  In  addition,  if 
we  consider  the  Lanchester  model  for  direct  fire  (the  square  law)  using  our  interpretation 
of  CB  and  CR  as  kill  rates  and  B{t)  and  R{t)  as  units  of  Blue  and  Red  value  at  time  t, 


i.e.: 


(3) 


dR{t) 

dt 


=  -CB  X  Bit) 


dBjt) 

dt 


=  -CR  X  Rit) 


then  r  =  y/CB  x  CRt  is  the  unit  of  non-dimensional  time  (see  Przemieniecki  [40]).  The 
quantity  y/CB  x  CR  reflects  the  rate  at  which  we  trace  the  attrition  path  -  a  clear  mea¬ 
sure  of  intensity.  (y/CB  X  CRt  appears  as  the  argument  in  the  cosh  and  sinh  terms  of  the 
solutions  of  (3).)  Finally  the  quantity  $o  =  is  known  as  the  superiority  param¬ 

eter  (see  Przemieniecki  [40]).  If  >  1?  Red  will  annihilate  Blue  before  being  annihilated, 
and  vice  versa  if  $o  <  1-  Thus  $o  is  the  initial  effective  force  ratio. 
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We  can  also  regard  the  potential/anti- potential  method  as  a  means  of  aggregation 
which  is  consistent  with  the  Lanchester  square  law.  Consider  the  Lanchester  system  of 
differential  equations: 


(4) 


dBm  dRi  dRn  \ 
dt  '  dt  '  '  dt  ) 


—  (5i,  .  .  .  ,  Bmt  Rli  •  ■  ■  1  dirij 


0 

L 


or 


/  dBi/dt 


0 


\dRn/dtJ  \Rn/ 

Here  we  have  assumed  that  at  all  times  t  the  kill  rates  are  constant.  If  we  fix  CB  and  CR, 
then  the  values  VBi  and  VRj  will  be  constant  up  to  some  common  positive  scale  factor. 
(Normalize  the  value  of  some  weapon  type,  say  VBi^  to  1  to  remove  this  ambiguity,  or 
carry  it  through  the  argument  below.) 

If  one  now  defines  the  aggregate  Blue  value  at  time  t  by 


m 

B{t)  =  Y^VBiX  B,{t) 

1=1 


and  the  aggregate  Red  value  at  time  t  by 


R(t)  =  J2VRj  xRj{t), 

j=i 

then  one  can  show  that  B{t)  and  R(t)  satisfy: 


(5) 


dt 

m 

dt 


=  -CR X  R 


=  -CB  X  B. 


The  proof  of  this  fact  is  straight-forward  (see  Section  Ic  below). 

This  result  says  that  given  the  eigenvalue  method  as  a  means  of  aggregation,  (5) 
is  a  consistent  aggregated  attrition  model  (see  Section  Ic  below).  In  effect,  aggregation 
commutes  with  attrition  —  we  can  run  the  attrition  model  (4)  for  a  period  of  time  and 
aggregate  the  result,  or  first  aggregate  and  then  run  the  attrition  model  (5)  -  either  way, 
the  result  will  be  the  same! 
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Note  that  the  converse  is  also  true.  If  we  seek  a  linear  aggregation 


B{t)  =  Y,VBiX  Bi{t) 

i=l 

n 

j=l 

where  Bi(t),  Rj{t)  are  the  solutions  to  (4)  and  VBi,  VRj  are  some  constant  “weights”, 
then  imposing  the  attrition  model  (5)  on  the  aggregated  forces  compels  the  relations 

n 

CB  xVBi  =  J2KijVRj 

j=l 

m 

CR  X  VRj  =  ^LjiVBi 

i=l 

i.e.  the  weights  must  be  the  ones  derived  from  the  eigenvalue  method.  (Again,  see  Sec¬ 
tion  Ic  below.) 

We  conclude  this  section  with  a  number  of  additional  observations: 

1)  li  Lx  K  and  K  x  L  are  irreducible  (as  required  by  the  Frobenius- Perron  theorem) 
then  one  can  show  that 


>  0  and 

uKii 


Thus  increasing  even  one  kill  rate  (holding  the  others  constant)  will  increase  Amax  snd 
our  measure  of  intensity  I  =  \/Xxna.x-  This  is  a  derived  result  which  lends  credence 
to  the  eigenvalue  method  and  our  interpretation  of  7  as  intensity.  For  numerical 
computation  of  these  derivatives  see  Section  4. 

2)  We  have  included  in  the  last  part  of  Section  1  a  geometric  proof  of  the  Frobenius- 
Perron  theorem  that  uses  the  Brouwer  Fixed  Point  Theorem.  This  theorem  says  that 

n 

every  continuous  map  from  the  closed  n-ball  =  {(a:i, . . .  ,a:„)  G  R”  ] 

i=l 

into  itself  has  a  fixed  point.  This  leads  to  the  consideration  of  more  general  fixed  point 
techniques  as  a  basis  for  assigning  force  values  and/or  aggregating  force  components. 
Fixed  point  methods  involving  contractive  maps  often  lead  to  iterative  procedures 
which  permit  the  effective  computation  of  the  fixed  point.  This  is  the  case  with  the 
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eigenvalue  method,  and  as  we  shall  see,  it  is  also  the  case  with  the  ATCAL  weapon 
importances. 

3)  When  the  Kij  and  Lji  depend  on  other  variables,  including  possibly  Bi  ,  .  .  .  ,  BjYlt 
Ri,...,Rn,  it  is  important  to  try  to  understand  the  resulting  dependence  of  V  Bi  and 
VRj  on  these  variables.  This  is  especially  interesting  in  a  stochastic  setting  where 
Kij,  Lji  are  random  variables.  We  can  ask  for  the  distribution  of  the  VBi,  VRj.  We 
have  investigated  these  questions  for  small  examples  and  the  results  are  presented  in 
Section  4. 

4)  The  potential/anti-potential  method  has  come  under  various  criticisms. 

a)  First,  it  is  clear  that  the  values  VBi,  VRj  are  really  measures  of  relative  lethality 
and  that  vulnerability  is  not  accounted  for.  That  is  because  they  are  instanta¬ 
neous  measures.  If  the  Kij,  Lji  depend  on  time,  perhaps  via  a  dependence  on 
the  current  force  structure  Bi,  Rj,  then  the  values  will  vary  over  time.  It  is  in 
this  change  over  time  where  vulnerability  enters. 

b)  “. . .  The  numeric  values  of  the  scores  are  sometimes  over  sensitive  to  small  changes 
in  the  input  kill  rate  matrices.  Zero  score  values  sometimes  occur  for  major 
weapons.  Also,  the  method  sometimes  splits  an  engagement  into  two  disconnected 
separated  engagements  . . .”. 

All  of  these  problems  stem  from  the  failure  of  LK  and  KL  to  be  irreducible  (or  from 
numerical  problems  that  are  the  result  of  being  close  to  reducible).  Simple  examples 
where  such  problems  occur  are  given  by; 


b.  Relationships  with  Lanchester  Attrition  Models  and  Aggregation. 

In  previous  reports,  we  have  shown  that  aggregation  cannot  be  considered  in  isolation 
from  the  attrition  model.  It  is  only  by  linking  the  two  that  we  can  account  properly  for 
the  dynamics  of  force-on-force  dependence  in  any  theory  of  combat  values.  To  illustrate 
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a  case  where  compatibility  can  be  achieved,  consider  the  potential/ anti-potential  method 
explained  above.  We  have 


( 


dBx 

dt 


\ 


dB_j^ 

dt 

dR, 

dt 


dRn 

dt 


0  -L^' 

-K'^  0 


Bm 

Rl 

\Rn/ 


aggregating  to 


—  =  -CRxR 
dt 


m 

dt 


=  -CB X  B 


where  B  =  Bi  and  i?  =  V Rj  x  Rj  is  the  aggregation  process. 

1=1  i=i 

If  we  assume  that  the  Kij  and  Lji  are  constant,  then  the  values  VBi  and  VRj  are 
constant  and  the  Jacobian  J  of  our  aggregation  process  is  the  2  by  m  -f-  n  matrix 

r_(VBi  ...  VBn  0  ...  0  \ 

0  ...  0  VRi  ...  VRnJ' 

Now  suppose  our  initial  force  composition  is  (5“ , . . . ,  B^,  i?2,...,il°)^and  that 

ip(t)  =  (Blit), . . . ,  Bm{t),  Ri{t), . . . ,  Rmit)f 
is  the  “path”  of  attrition  in  the  state-space  R"‘+”.  This  means  that  (f{t)  satisfies 

0  -L^' 

-K'^  0 


dip 

dt 


and 


lO  \r 


V(0)  =  (B»,...,O 


(These  conditions  uniquely  determine  (/?.)  Now  consider  the  projection  of  this  path 


into  the  state-space  R^.  It  satisfies 


dj 

dt 


dBi 


T 


(This  is  where  the  fact  that  Kij  and  Lji  are  constant  comes  into  play.)  In  other  words, 

c?7  d(^ 

dt  dt 

which  is  the  standard  relationship  between  the  velocity  vectors  of  a  path  and  those  of  its 
image  under  some  map.  On  the  other  hand, 

Id" 

\-CB  0  / 

because  of  our  choice  of  attrition  model  for  the  aggregated  forces.  The  compatibility  we 
have  been  discussing  is  the  requirement  that 

(  ° 

\  -CB  0  /  ^ 

For  example,  we  must  show 

^  dBi 

-CR  X  i?  =  ^ 

i=l 

But  the  right-hand  side  is  ^  -  just  differentiate  the  aggregation  equation 
B  =  ^  BiXVBi.  Thus  the  projected  path  7(t)  is  the  path  of  attrition  in  the  aggregated 
variables  B  and  i2!  This  establishes  consistency  of  the  potential  anti-potential  method 
with  the  Lanchester  attrition  models, 
c.  The  Frobenius-Perron  Theorem. 

In  this  part  of  Section  1,  we  discuss  the  Frobenius-Perron  Theorem  which  plays  a  key 
role  in  the  development  of  the  so-called  potential/anti-potential  or  eigenvalue  method  for 
assigning  values  (scores)  to  weapon  types.  It  turns  out  that  a  very  similar  method,  based  on 
the  same  Brouwer  fixed  point  theorem,  can  be  used  to  establish  existence  and  uniqueness 
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for  ATCAL’s  weapon  importances.  We  conclude  our  discussion  with  some  observations  on 
computational  methods  which  implement  the  theorem. 

Let  A  be  a  real  n  x  n  matrix  with  aij  >  0  for  every  i,j: 


ail 

Ol2 

*  ^In  ^ 

021 

022 

•  ^2n 

V  Ujii 

On2  • 

•  *  ^nn  ' 

Such  an  A  is  said  to  be  non-negative. 

Definition  1.  A  matrix  A  is  said  to  be  reducible  if  the  index  set  {1,2,  ...,n}  can  be 
partitioned  into  two  disjoint  sets  {*i ,  •  •  ■ ,  v)  and  {ki, . . . ,  k^}  (p  -\-  u  =  n)  such  that 

^ia  kg  ~  0  Oi  1,...,^  fd  I,...,!/. 


In  other  words,  after  applying  a  suitable  permutation  to  both  the  rows  and  columns  of  A, 
it  takes  the  form 

p  u 


If  a  matrix  is  not  reducible  it  is 


V 
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said  to  be  irreducible. 


□ 


In  1907  Perron  proved  a  remarkable  theorem  concerning  the  eigenvalues  and  eigen¬ 
vectors  of  a  positive  matrix  (ojj  >  0).  This  was  later  generalized  by  Frobenius  to  the  case 
of  irreducible  non-negative  matrices. 

Theorem  2.  (Perron-Frobenius.)  An  irreducible  non-negative  matrix  A  =  (aij)  has  a 
unique  largest  real  positive  eigenvalue  Amax  which  is  a  simple  root  of  the  characteristic 

polynomial  x”  —  trA  x”~^  -1 - 1-  (-1)"  det  A  of  A.  The  modulus  (complex  absolute  value) 

of  any  other  eigenvalue  will  be  less  than  or  equal  to  Amax;  and  it  will  always  be  strictly  less 
than  Amax  is  positive.  To  this  maximal  eigenvalue  Amax  there  corresponds  a  positive 
eigenvector  v  =  (ui, . . . ,  with  u,-  >  0  for  every  i.  Of  course,  v  is  determined  only  up 
to  a  positive  scalar. 
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We  will  not  prove  this  result  in  its  full  generality  here.  The  reader  should  consult 
Gantmacher  [20]  for  complete  details.  Instead,  we  will  give  a  topological  proof  that  makes 
use  of  Brouwer’s  fixed  point  theorem. 

Let  .4  be  a  real  non-negative  irreducible  n  xn  matrix  and  let  x  =  (xj , . . . ,  Xn)^) 
Xi  >  0  all  i,  ^  ^  0  be  a  non-zero  non-negative  (column)  vector.  Since  the  theorem  is 
trivial  for  n  =  1,  we  will  assume  n  >2. 

Lemma  3.  4.x  7^  0. 

Proof:  7  has  at  least  one  strictly  positive  component,  say  Xj  >  0.  If  4x  were  to  be  0, 
we  would  be  forced  to  have  a,j  =  0  for  every  j  =  1, . . .  ,n,  but  a  column  of  zeros  would 
make  A  reducible,  contrary  to  our  assumption.  D 


Recall  that  the  standard  (n  -  l)-simplex  C  R”  is  defined  by 

=  |(xi,. . .  ,x„)'^  e  R"  s.t.  Xi  >  0  for  every  i  and  ^  x^  =  1 
For  example,  when  n  =  2,  AMs  a  line  segment  in  the  plane; 


t=i 


and  when  n  =  3,  A^  is  a  triangle  in  space: 


Proof  of  Theorem  2.  We  define  a  continuous  map  g:  A”  ^  — »•  A”  ^  by 


9i^) 


Ax 

a{Ax) 
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where  (T((yi, .  ..,ynY)=  ^Vi,  i-e.  a  means  add  the  entries  in  the  n- vector  to  which  it  is 
applied. 

Note  that  a:  R"  -»  R  is  continuous,  so  that  its  reciprocal  will  be  continuous  away 
from  places  where  it  is  zero.  In  fact,  Lemma  3  shows  that  a(Ax)  >  0  for  x  G  A”  Also 
note  that  g(x)  is  again  in  A””^  since  sr(i)  is  a  non-zero  non-negative  vector  with 


a{Ax) 

a{Ax) 


Now  Brouwer’s  Fixed  Point  Theorem  says  the  following: 

.  n 

Let  B”  =  {x  €  R"  s.t.  |xp  =  ^  xf  <  1}  be  the  unit  n-ball  in  R"  consisting  of  all 

i=l 

vectors  of  length  less  than  or  equal  to  1  in  R",  and  let  /:  B"  B"  be  any  continuous 
map,  then  /  has  a  fixed  point,  i.e.  there  exists  a  vector  x  G  B”  such  that  /(x )  =  x. 


This  theorem  applies  as  well  to  any  topological  space  which  is  topologically  equivalent 
(homeomorphic)  to  B”.  This  includes  A”.  Thus  our  map  g:  A”  ^  — >  A"  ^  has  a  fixed 
point  X  G  A"“^.  We  have 

Ax  — 

- =r-  =  X, 

cr{Ax) 


or,  setting  A  =  cr(A.x)  >  0, 


Ax  =  Ax. 


In  other  words  A  >  0  is  an  eigenvalue  for  A  with  non- negative  eigenvector  x .  In  fact,  every 
component  of  x  must  be  strictly  positive.  To  see  this,  suppose  to  the  contrary  that  x  has 
some  0  entries.  Without  loss  of  generality,  we  can  permute  the  coordinate  axes  so  that 
^  =  (xi , . . . ,  Xfc,  0, . . . ,  0)^  with  Xi  0  for  i  =  1, . . .  ,k  where  1  <  k  <  n.  We  then  have 


However,  this  relation,  plus  the  fact  that  x\  >  0, 
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X2  >  0,...,Xk  >0,  forces  A21  to  be  an  {n  -  k)  x  k  block  of  zeros,  making  A  a  reducible 
matrix. 

We  have  thus  established  that  a  real  non-zero  non-negative  irreducible  matrix  A  has 
a  positive  eigenvalue  with  corresponding  positive  eigenvector.  Additional  arguments  are 
needed  to  show  that  this  eigenvalue  has  the  largest  absolute  value  among  the  real  eigen¬ 
values  of  A,  and  that  no  other  real  eigenvalue  has  a  positive  eigenvector.  Notice  that  the 
fixed  points  of  g  are  in  one-to-one  correspondence  with  the  “non-negative”  eigenvectors 
"y  =  (yi , . . . ,  y„)^,  yi  >  0  all  i,  of  A  (up  to  positive  scale),  and  that  the  argument  above 
shows  that  all  the  fixed  points  lie  in  the  interior  of  A"''\  which  is  equivalent  to  yi  >  0  for 
all  i.  The  theorem  of  course  implies  that  g  has  a  unique  fixed  point. 

We  remark  that  reducible  non-negative  matrices  will  still  have  non-negative  eigenval¬ 
ues  with  non-negative  eigenvectors,  however  there  may  be  several  such.  Examples  to  think 
about  are: 


1 

0 

1' 

'1  O' 

1  O' 

[0 

0_ 

7 

0  1 

7 

0  2 

Computational  Issues 

One  can  compute  the  unique  fixed  point  of  the  map  g:  A”“^,  where  y(x)  = 

,  by  an  iterative  procedure.  For  this  to  work,  we  must  assume  that  A,  in  addition  to 

(t{Ax  ) 

being  non-negative  and  irreducible,  is  primitive.  If  all  the  entries  of  A  are  positive,  then  A 
is  primitive.  If  A  has  some  zero  entries,  then  A  will  be  primitive  if  and  only  if  some  power 
of  A  is  positive. 

To  compute  the  fixed  point,  we  pick  any  initial  vector  6  A”~\  say  = 

. . . ,  Xn  where  xj^^  >  0  all  j  =  1, . . . ,  n  and  ^  =  1,  and  set 

=  y(x^*^) 


so  that 

x^^^  =  x^^^  =  g(x^^^)  =  g(g(x^^^)) 

and  so  on.  One  can  show  that: 
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Proposition  4.  lim  x^''>  =  xo  where  g{xo)  =  xq  is  the  unique  fixed  point  of  g  and 

i—^oo 

therefore,  up  to  scale,  the  unique  positive  eigenvector  of  A. 


This  iterative  method  breaks  down  if  A  is  reducible  or  if  A  is  not  primitive,  as  the 
example: 

shows. 


Sensitivity 

Let  A  be  a  non-negative  irreducible  matrix.  The  largest  eigenvalue  Amax  can  be 
characterized  as  follows: 


A 


max 


max 


(min 

l<i<n 

Xi^O 


where  the  maximum  is  taken  over  all  vectors  x  =  (xi , . . . ,  with  Xi  >  0  for  all  i 

n  ^  ^ 

and  ^  Xi  =  1  and  where  (Ax)i  =  aijXj  is  the  component  of  the  vector  Ax  (see 
j=i  >=1  ^ 

Gantmacher  [  ]).  Moreover,  the  unique  vector  x^°^  in  giving  the  maximum  is  in  the 

interior  of  A”~^  and  is  an  eigenvector  with  eigenvalue  Amax: 


Ax^  ^  —  Ajnax  ^ 


(0) 


Now  suppose  every  element  in  the  column  of  A  is  increased  by  £  >  0  yielding  a 
new  non-negative  irreducible  matrix  Aj^s-  We  compute  that  (Aj,gX^°^)i  =  (Ax^*^^)i  -|-£Xj 
so  that 

{Aj,eX^^'>)i 


mm 

l<i<n 


>  Aniax  +  min  £ 


(0) 
T  • 

l<t<n  3;^®^ 


It  follows  that  Aniaxi  the  largest  eigenvalue  for  Aj^g  is  at  least  A^ax  +  £  ^  ^max- 

Now  suppose  in  addition  that  A  is  primitive.  From  the  definition  of  primitive,  it 
follows  that  A^  =  {a\f)  has  only  positive  entries  for  some  q  >  1.  If  we  increase  the  i,j- 
entry  of  A  by  £,  then  every  entry  in  the  j*'*'  column  of  =  A’ A  increases  by  at  least 

(  min  a^ki)s.  The  argument  above  then  shows  that  the  largest  eigenvalue  of  A^'^^  (which 

l<k<n 
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is  A^+^)  increases  by  at  least  cqC  where  co  is  fixed  positive  constant  independent  of  £.  It 
follows  that 


Co  < 


=  (g  +  IjAJ, 


doi 


dai 


tj 


so  that  is  positive.  Moreover,  we  can  give  a  positive  lower  bound  for  this  value. 

oaij 

Numerical  confirmation  of  this  result  can  be  found  in  Section  4. 
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§2.  Fire  Allocation  Methods  for  Assigning  Weapon  Values 
a.  The  Methods. 

i.  Introduction  to  Fire  Allocation  Methods. 

Fire  allocation  methods  for  assigning  values  to  the  weapons  present  in  an  engagement 
are  based  on  the  principle  that  the  value  of  opposing  weapons  at  an  instant  in  time  will 
determine  (along  with  other  factors/variables)  how  we  allocate  the  fire  of  friendly  weapons. 
This  allocation  of  fire  simultaneously  impacts  the  rate  of  attrition  of  various  weapon  types, 
so  we  will  need  attrition  models  that  capture  this  effect.  Finally  attrition  rates  (at  a  given 
instant  in  time)  effect  value  at  that  instant,  because  we  will  assume,  as  we  did  in  the 
eigenvalue  method  (see  §1),  that  the  value  of  a  weapon  is  proportional  to  the.  rate  at  which 
it  is  destroying  enemy  value. 

As  one  might  expect,  working  with  models  and  equations  that  portray  such  a  com¬ 
plex  interaction  between  value,  fire  allocation,  and  attrition  can  be  difficult.  Unlike  the 
relatively  simple  linear  equations  and  matrix  techniques  of  the  eigenvalue  method,  fire  al¬ 
location  methods  invariably  lead  to  non-linear  equations  with  parameters.  In  this  section, 
we  will  illustrate  fire  allocation  methods  with  examples  that  use  a  Lanchester-like  attrition 
model.  Other  attrition  models  could  be  used.  In  fact,  in  a  very  general  sense,  the  ATCAL 
model  (see  §3)  also  embodies  a  fire  allocation  method. 

Because  the  examples  lead  to  complex  systems  of  non-linear  polynomial  equations, 
we  also  spend  some  time  in  this  section  discussing  methods  to  solve  these  equations.  For 
relatively  few  weapon  types  on  each  side,  such  systems  of  non-linear  equations  can  be 
solved  with  the  aid  of  a  mathematical  computation/symbolic  computation  package  such 
as  Mathematica  or  Maple.  Finally,  we  can  observe  that  despite  the  presence  of  non-linear 
equations  in  our  fire  allocation  examples,  the  method  exhibits  features  that  are  surprising 
similar  to  the  eigenvalue  method^ 

ii.  A  Simple  Model. 

We  first  consider  a  simple  example  with  two  weapon  types  on  each  side.  We  denote 
by  2>i(t),  f>2(f)  the  number  of  blue  weapons  of  type  1  and  type  2  respectively,  present  in 
the  engagement  at  time  t.  Likewise  ri(t),  r2{t)  are  the  respective  numbers  of  each  type  of 
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red  weapon  present  at  time  t. 

Our  attrition  model  will  be  a  simple  Lanchester  square  law  based  on  the  ordinary 
differential  equations: 

dhi  ,  . 

—  =  -pi(ciiri  +  Ci2r2) 

— ^  =  — P2(C21?'1  +  ^22^2) 
at 

^  =  -qi(anh  +  012^2) 
at 

=  —52(021^1  +  0.22^2)- 

The  Cij  and  aij  are  attrition  rates  that  reflect  the  rate  at  which  one  unit  of  a  particular 
weapon  destroys  units  of  another.  Thep*;  and  represent  an  allocation  of  fire.  Specifically, 

we  will  have  pi,P2,  9i,?2  >  0  with  pi  4-p2  =  1  and  qi  +  q2  =  1- 

In  this  particular  model  we  are  assuming  both  red  weapon  types  fire  at  a  given  blue 
target  in  the  same  proportion  and  likewise  for  blue.  A  more  realistic  model  would  perhaps 


—  =  -CiiPiiTi  -  Ci2P12J’2 

at 

db2 

— —  =  — C2iP21»'1  —  C22P22f'2 
at 

-ry  =  — ~  <212912^2 
at 

— ^  =  —021921^1  ~  <*22922^2 
at 

where  pkf  >  0,  qkt  >  0  with  pik  +  P2k  —  1  for  A;  =  1, 2  and  qn  +  q2t  =  1  for  ^  =  1,2. 

In  general  pi,p2,9i,92  will  be  functions  of  various  parameters.  To  illustrate  the 
method,  we  select  an  allocation  of  fire  that  assumes  we  fire  at  each  target  in  propor¬ 
tion  to  its  value.  Thus  if  a  unit  of  61  is  twice  as  valuable  as  a  unit  of  then  2/3  of  red’s 
fire  will  be  directed  at  type  hi  and  1/3  at  type  62-  Specifically,  we  set 

hM^) 

hibi{0)  -f-  /i2^2(0) 

_ ^2^2(0) 

^  "b  ^2^2(0) 

91^1  (0) -h  92^2(0) 

_ 92^2(0) 

9iri(0) -b  92^2(0) 
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where  hi  is  the  “value”  of  one  unit  of  6i  at  the  current  time  (which  we  take  as  <  =  0  for 
simplicity)  and  /12  is  the  “value”  of  one  unit  of  62.  Thus  /ii6i(0)  +  h2b2{0)  is  the  total 
“value”  of  blue  at  the  current  time.  Similarly  gi ,  52  ^re  the  values  for  red  at  the  current 
time. 

This  may  or  may  not  be  a  “realistic”  model  of  fire  allocation.  What  is  important  is 
that  the  fire  allocation  factors  {pk  and  qe  in  this  case  are  functions  of  the  unknown  values 
and  other  known  quantities  (fct(O)  and  rj(0)  in  this  case).  For  a  discussion  of  the  general 
case,  see  Section  2.a.v  below. 

As  with  the  eigenvalue  method,  we  assume  that  value  is  proportional  to  the  rate  at 
which  a  weapon  destroys  opposing  value.  For  example, 

i?iri(0)  \  ,  f  ^2^2(0) 


Xhi  -  gi 


an  +  92 


fifiri(O)  +  g2J’2(0) 


021 


,5l''l(0)  +  i/2»'2(0), 

which  says  that  the  value  of  one  unit  of  61  is  proportional  to  the  rate  at  which  it  destroys 
red  value.  Note  that  the  units  on  the  constant  of  proportionality  A  are  (time)  .  This 
equation  can  be  viewed  as  being  derived  from 


Xhibi  —  gi 


9iri{0) 


firiri(O)  +  j?2r2(0) 


anbi  +  fl'2 


3^2  ^2(0) 


021^1 


.fi^in(O)  -  fir2r2(0) 

which  says  that  the  total  value  of  all  weapons  of  type  bi  is  proportional  to  the  rate  at 
which  they  collectively  destroy  red  value  at  the  current  time. 

This  leads  to  the  following  system  of  non-linear  equations: 

qiigi^i(0)  + 021^^2^2(0) 

5in(0)  +fir2»'2(0) 

012^1^(0)  +  022^2^2(0) 
giri(O)  -f-  ^2^2(0) 

Cii/li6i(0)  -f  C2i/l2^2(0) 

^i^i(O)  +  ^*2^2(0) 
ci2^jfei(0)  +  ^22^2^2(0) 

^i^i(O)  +  ^2^2(0) 


Xhi  = 
Xh2  = 
Xgi  = 
Xg2  = 


or 


Xri{0)higi  -I-  Xr2(0)hig2  -  aiiri(0)5i  -  a2ir2(0)5|  —  0 
Ari(0)/i2fifi  +  Ar2(0)/i25'2  ~  Oi2^i(0)fifi  ~  022^2(0)^2  “  0 
Xbi{0)higi  +  Xb2(0)h2gi  —  Cii6i(0)/ii  —  02162(0)^2  =  0 
A61  (0)61^2  +  Xb2{0)h2g2  —  01261(0)61  —  02262(0)62  =  0. 
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We  see  that  for  each  A,  we  have  a  system  of  four  homogeneous  quadratic  equations  in 
four  variables  9i,g2-  In  general  we  expect  such  a  system  to  have  no  solutions 

other  than  the  trivial  one  hi  =  0,  hi  =  0,  gi  =  0,  $-2  =  0-  However  for  certain  values 
of  A  (the  analogue  of  an  eigenvalue  in  the  linear  case)  there  will  be  non-trivial  solutions. 
The  problem  is  to  determine  what  values  of  A  >  0  cause  the  associated  system  to  have 
a  solution  with  hi  >  0,  hi  >  0,  gi  >  0,  92  >  0.  Thus  our  fire  allocation  method  is  a 
non-linear  analogue  of  the  eigenvalue  method.  Note  that  value  is  determined  only  up  to 
scale,  so  that  we  must  normalize  the  value  relative  to  one  type  (which  then  has  value  1). 
There  are  two  procedures  for  solving  this  problem.  The  first  treats  the  whole  system  of 
equations  and  uses  classical  elimination  theory  to  compute  a  polynomial  iZ(A)  whose  roots 
are  the  values  of  A  for  which  the  system  has  a  non-trivial  solution.  The  coefficients  of  i2(A) 
are  polynomial  expressions  in  a^j,  Cji,  6^(0)  and  rj(0).  The  specific  computational  device 
for  finding  i2(A)  is  known  as  the  Macaulay  resultant  and  it  is  discussed  below. 

A  second  approach  leads  to  two  sets  of  equations  analogous  to  a  simultaneous  eigen¬ 
value  problem 


Pi(A,  hi,  hi)  =  0 

P2(A,/li,/l2)  =  0 


Qi(A,yi,^2)  =  0 

Q2i^i  91,92)  =  0 


where  the  first  set  of  equations  Pj  =  0,  P2  =  0  involves  only  blue  values  and  A  and  the 
second  set  Qi,  Qi  involves  only  red  values  and  A.  These  equations  are  quartic,  homogeneous 
in  hi ,  hi  or  gi ,  92  and  are  derived  by  substitution.  For  example,  using 


aiiri{0)gi  +  021^2(0)^2 
^iri{0)gi  +  Xri{0)g2 

qi2ri(0)gi  -b  022^2(0)^2 

Ari(0)gi  +  Xri{0)gi 
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and  substituting  into  the  equations  for  gi  and  g2  yields: 

Affi  = 

\g2  = 

Clearing  denominators  gives 

0  =  A^g'i(ri(0)sfi  +  r2{0)g2)  [(aii6i(0)ri (0)  +  012^2(0)^1(0))^! 

+(o21^l(0)j’2(0)  +  O22^2(0)?’2(0))5r2] 

-  Cii6i(0)(aiiri(0)firi  +  a2ir2(0)5r|)^  -  C2i62(0)(ai2»'i(0)fir?  +  0,22^2(0) gjf 
0  =  g2i'i'i{0)gi  +  r2(0)fir2)  [(aiifei(O)ri(O)  +  ai2i2(0)ri(0))5ri 
+  (021^1(0)^2(0)  +  022^2(0)?'2(0))5'2] 

—  Ci26i(0)(aiiri(0)s'i  +  O2ir2(0)fir|)^  —  C22^»2(0)(oi2ri  (0)5^1  +  O22^2(0)5r|)^. 

The  right  had  sides  are  Qi{^,gi,g2)  and  Q2{^y9i^  92)  respectively. 

Again,  the  system 

Pl(A,/li,/l2) 

P2(A,  hi,h2) 

generally  has  no  solutions.  However  for  certain  values  of  A  it  will.  Likewise  for 

Qi(A, 91,52)  =  0 

Q2(A,9i,92)  =  0. 

In  each  case  there  is  a  polynomial  in  A,  Rp(\)  or  Pq(A),  whose  roots  are  those  A  for  which 
that  particular  system  has  a  non-trivial  solution.  Optimally,  we  will  find  a  unique  positive 
A,  call  it  Aj„3_x,  ivith  \ 

■^p(Amax)  ~  0  and  PQ(Ainax)  —  0 

which  will  lead  to  a  unique  solution  for  the  weapon  values  (after  normalizing).  We  cannot 
as  yet  prove  this. 


>  /'n\  t  aiin(0)5?+a2ir2(0)gn 
Cll&ljUJ  ^  Ari(0)gi  +  Ar2(0)g2  ) 

1  { ail>'l(0)g?  +  °2i  *~2(0)g2 

Ari(0)gx  +  Ar2(0)g2 


u  ic\\  ( ai2n(o)g?+a22''2(o)g2 
+  C2102(0)(^  Ari(0)giH-Ar2(0)g2 

A  ,  L  /r>\  /ai2*'i(0)g?+922r2(0)g^A 
)  +  ^2(0)  [  Arx(0)g!  +  Ar2(0)g2  j 

-  7  / a  12 »~i(0)g?  + 022^2(0)32 

I  +C2202tUjl  Ari(0)gi  +  Ar2(0)g2 


Cl2il(0)(2^ 


iin(o)g? +a2i^2(Q)g2  ^ 
Ari(0)^i  +  Ar2(0)p2  J 

aiiri{0)g^-\-a2ir2{0)9l 
Ari(0)yi  +  Ar2(0)^2 


^  +^22^2(0)^^ 


+  62(0) 


ai2n(Q)g?  +  g22r2(Q)g2 

Ari(0)^i  +  Ar2(0)52 
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iii.  An  Iterative  Scheme  with  Shifting  Fires  and  Changing  Values. 

Before  giving  a  specific  numerical  example  and  discussing  more  general  fire  allocation 
schemes  for  determining  weapon  values,  we  note  that  the  simple  method  above  leads  to  an 
interesting  iterative  time-step  scheme  that  links  value,  attrition,  and  fire  allocation. 

At  time  to  with  forces  6i(to),  &2(<o),  ri{to),  r2(to)  present  in  the  engagement,  we 
calculate  the  current  values  hi(to),  h^ito)-,  gi{to):  5'2(^o)-  With  these  values  known  we  can 
update  the  attrition  for  a  time  step  At  by  setting 

^i(^o  +  =  bi(to)  -f-  A6i 


etc.,  where  the  attrition  over  the  time  step  is: 


A6i 


hi{tQ)bi{to) 


hi{to)bi{to)  +  h2{to)b2{to) 


(aiiri(to)  +  ai2r2(to))^^ 


etc.  This  updates  the  forces  present  at  time  t\  —  to  T  We  now  cycle  through  the 
process  again.  With  this  sort  of  model,  we  expect  to  see  values  change  and  fires  reallocate 
as  losses  occur. 

To  obtain  a  time  independent  value  for  a  unit  of  say  6i,  we  could  average  its  value 
over  time 

1  r^final 

hi  =  - — —  /  hx{t)dt 

*final  ^init 


or 


nAi  +  At)  -I - 1-  /li(tinit  +  (w  —  l)At)  -f-  2^l(*init  +  nAt'^ 

where  nAt  =  tfinal  ~  ^init  is  the  duration  of  the  engagement  which  begins  at  time  tjnit  and 
ends  at  time  tfinal-  Alternatively,  we  could  compute  a  somewhat  different  average 


hi  = 


r  tfinal 
*'  ^init 


bi{t)hi{t)dt 


which  measures  the  average  contribution  to  blue’s  total  value  by  a  single  unit  of  b^  at 
the  start  of  the  engagement  over  the  duration  of  the  engagement.  Weapons  with  high 
vulnerability  will  not  fare  as  well  as  weapons  of  equal  lethality  but  less  vulnerability  with 
this  measure  of  value. 
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iv.  A  Numerical  Example  of  the  Procedure  for  Solving  the  Simultaneous  Sys¬ 
tem  of  Non-linear  Equations  Arising  in  our  Simple  Fire  Allocation  Model. 

Let’s  take  our  initial  force  strength  to  be 

6i(0)  =  100  ri(0)  =  150 

62(0)  =  40  r2(0)  =  50 


with  attrition  coefficients 

ail  =  -05  ai2  =  -03 

021  —  -02  022  — 

cii  =  .06  C12  =  .04 

C21  =  .03  C22  =  -08 

so  that  the  attrition  model  becomes 


dt 

dt 

dri 

dt 

dr2 

dt 


/  lOO/ii  \ 

VlOOhi  +40/12  J 
/  40/t2  \ 

Vioo/ii  +40/12/ 

(  150gi  \ 

V  150^1  +  50g'2  / 

(  50ff2  \ 

VlbOjfi  +  50gf2/ 


(.06ri  +  .04r2) 
(.03ri  +  08. r2) 
(.0561  +  .0362) 
(.0261  +  .0762). 


Note  that  the  red  weapons  have  the  better  kill  rates  and  can  be  expected  to  have  better 
values.  Thus  gi ,  g2  should  be  larger  than  hi ,  /12  respectively. 


Our  equations  for  the  values  are 


1) 

2) 

3) 

4) 


^  7.5gf  +  l.OgI 
ISOgfi  +  50g'2 
_  4.5^1  +  3.5^2 
150^1  +  505^2 

6.0/ij  +  1.2/i2 

^  lOO/ii  +  40/12 
AMl  +  3.2hl 
~  loo/ii  +  40/12  ■ 


After  clearing  denominators,  we  are  led  to  a  system  of  four  homogeneous  quadratic  equa¬ 
tions  in  four  variables  hi,h2,  51,5^2-  We  treat  A  as  a  coefficient  rather  than  a  variable.  The 
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Macaulay  resultant  of  this  system  will  be  a  polynomial  in  A  whose  roots  are  those  values 
of  A  for  which  the  system  has  a  solution.  (See  Section  2b  for  details  on  the  Macaulay 
resultant.) 

An  alternative  approach  that  is  computationally  more  tractable,  is  to  solve  equations 
1)  and  2)  for  hi  and  /i2  and  substitute  those  expressions  into  3)  and  4).  Likewise  we  can 
solve  3)  and  4)  for  gi  and  and  substitute  into  1)  and  2).  This  leads,  after  some  simple 
algebra,  to  two  systems  of  two  equations  each  linked  by  A.  Specifically,  if  we  let  u  =  hi/ h2 
and  V  =  gi! g2i  'ffQ  arrive  at 

5)  0  =  (-286  +  lOOOA^)^^  +  (ddOOOA^)^^  +  (-133.6  +  34000A2)u2 

+  (13600A2)u  + (-21.04) 

6)  0  =  (-218)u'‘  +  (IIOOOOA^)^^  +  (-154.4  +  44000A2)u2  +  (34000A2)u 

+  (-42.32+  13600A2) 

and 

7)  0  =  (-361.8  +  139500A2)t;^  +  (45600A2)t;^  +  (-127.8  +  36000A2)t;2 

+  (12000A2)u  +  (-20.7) 

8)  0  =  (-2S9.S)v*  +  (139500A2)v3  +  (-160.8  +  dOSOOA^)^^  +  (36000A2)?; 

+  (-43.2  +  12000A2). 

Note  that  (5)-(8)  is  very  similar  to  the  simultaneous  eigenvalue  problem  that  arises  in  the 
potential/antipotential  or  eigenvalue  method  in  Section  1. 

Now  consider  equations  (5)  and  (6).  For  a  fixed  value  of  A  we  have  two  polynomials 
in  the  single  variable  u.  In  general  they  will  not  have  a  root  in  common,  but  for  certain 
A’s  they  will.  These  values  are  found  by  computing  the  resultant  (an  8x8  determinant  in 
this  case)  of  the  two  polynomials  (see  Section  2b).  This  resultant  (up  to  a  constant  factor) 
is: 

0  =  4.13749  -  16328A^  +  2.66084  x  10^ A^  -  2.22584  x  10^° A®  +  9.42028  x  10^^  A® 

-  1.5669  X  10^^A^°. 

The  only  positive  root  is  A  =  .0467053. 
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Similarly  taking  the  resultant  of  (7)  and  (8)  we  get  another  polynomial  in  A  which 
also  has  A  =  .0467053  as  a  root!  Taking  this  common  value  of  A,  we  find  the  common  root 
of  (5)  and  (6)  to  be 

M  =  /ii//i2  =  1.09535 

and  the  common  root  of  (7)  and  (8)  to  be 

V  =  Qx!  g2  —  1.05002. 


Normalizing  /i2  to  1,  we  get 

hx  =  1.09535 

hi  =  1.00000 

gx  =  1.20256 
gi  =  1.14534. 

Note  that  there  is  error  in  these  computations  and  that  the  results  are  accurate  only  to 
five  significant  figures. 


V.  General  Fire  Allocation  Methods. 

The  simple  example  presented  in  section  2.a.ii  above  can  be  easily  generalized.  For 
example,  we  could  consider  a  more  general  Lanchester  attrition  model 


=  -'^PikCikrkit)  -  ^p'iec'iibi{t)re{t) 
k=i  e=i 


r=l 


3=1 


i  =  1, . . . ,  m 
;■  =  l,...,n 


that  includes  both  the  square  and  linear  laws  (direct  and  indirect  fire).  The  fire  allocation 

factors  pik,  p'if,  qjri  and  9'  ^  must  satisfy 

\ 

n  n 

+  =  l  i  =  l,...,m 

j=i  j=i 


and 

m  m 

+  1  i  =  l,...,n. 

t=i  t=i 
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These  factors  can  be  functions  of  numerous  parameters  including  the  kill  rates  aji,  aji,  Cij, 
and  c[j,  the  force  levels  bi(to),  s-t  the  current  time  t  =  to,  and  other  constant  factors 

(for  example  the  analogues  of  the  availability  factors  in  ATCAL).  Most  importantly  they 
should  involve  the  unknown  weapon  values,  hi  and  gj,  of  the  blue  weapons  of  type  i  and 
the  red  weapons  of  type  j,  respectively.  These  values  hi,gj^  once  determined,  should  be 
interpreted  as  values  at  the  current  time  hiito),  gj{io)- 

We  are  assuming  here  that  each  weapon  type  mounts  one  weapon.  If  this  is  not  the 
case,  then  that  particular  type  must  appear  more  than  once  on  the  right  hand  side  of  our 
attrition  equations.  Moreover,  if  a  particular  weapon  is  not  capable  of  both  direct  and 
indirect  fire,  it  would  not  appear  in  both  places,  and  the  corresponding  allocation  factors 
would  not  appear.  For  example,  the  system  below: 


dt 

dh 

dt 
dr  I 
dt 
dr2 
dt 


=  -Piiaiin  -Pna'iihri  -  p'i2a'i2^i^2 

=  — P21<^21^1  —  P21*^21^2^1  ~  P22^22^2^2 
=  —  911^11^1  +  9i2^12^1^2 

=  —521^21^1  “  921^21^1  +  $22^22^2^2 


represents  a  red  type,  rj ,  that  mounts  a  weapon  capable  of  both  direct  and  indirect  fire, 
a  red  type,  r2,  that  mounts  an  indirect  fire  weapon,  a  blue  type,  fei,  that  mounts  two 
separate  direct  fire  weapons,  and  a  blue  type,  62?  that  is  an  indirect  fire  weapon. 

As  always,  we  assume  that  a  weapon’s  value  is  proportional  to  the  rate  at  which  it 
destroys  opposing  value;  so  that 


rt  n 

^h\  =  3j‘''j  i^o ) 

i=i 


n  n 

Xhjfi  =  ^  QjmO'jmQj  ^ 
j=l  j=l 
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m  m 

Ayi  =  +  J^PiiCiihibi(io) 

1=1  1  =  1 

m  m 

1=1  1=1 

This  is  a  system  oi  m  +  n  equations  in  m  +  n  variables  /ii,  .  .  .  ,  hm,'}  gi,...,gni  treating  A 
as  a  parameter  in  the  coefficients.  It  is  non-linear  in  general,  because  the  expressions  for 
qjr,  q'js,  Pik,  and  p'ii  will  involve  hi, . . .  ,hm,  gi,  ■  ■  ■  ,gn-  Moreover  the  equations  should  be 
homogeneous,  since  scaling  all  the  values  by  the  same  amount  should  not  cause  a  change 
in  the  allocation  of  fire,  because  the  relative  values  will  remain  the  same.  To  solve  systems 
of  this  sort  see  2b  below.  Solving  the  system  in  this  case  means  finding  those  A  for  which 
the  system  has  a  non-trivial  solution.  Again  this  is  done  with  resultants. 

Fire  allocation  methods  are  not  limited  to  use  with  Lanchester  attrition  models.  Any 
attrition  model  is  a  candidate  for  this  method  as  long  as  we  can  link  attrition  rates,  values, 
and  allocation  of  fires  as  we  have  done  above. 

b.  Resultants  and  Non-linear  Equation  Solving  Techniques 

This  section  is  meant  to  serve  as  an  introduction  to  resultants.  Simple  examples 
are  presented  which  illustrate  each  technique.  The  reader  should  have  little  difficulty  in 
applying  the  various  methods  to  any  particular  problem. 

Resultants  are  used  to  solve  systems  of  non-linear  polynomial  equations,  to  determine 
whether  or  not  solutions  exist,  or  to  reduce  a  given  system  to  one  with  fewer  variables 
and/or  fewer  equations. 

The  typical  input  will  be  a  system  of  m  equations  in  n-variables: 

)  •  •  •  r  ^n) 

f m(^^l  ^  ^  ^n)  0* 

Each  equation  has  an  associated  degree  d,  >  1.  Recall  that  /t(xi , . . . ,  has  degree  di 

n 

if  all  monomials  appearing  in  fi  have  X]  least  one  monomial 

1=1 
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has  ^  Ci  =  di.  Example:  /(xi,X2,X3)  =  Zx\xz  +  4xiX2  —  X2  +  7xz  —  1  has  degree  d  —  Z. 
The'integers  m,n,di, . . .  ,dm  are  important  indicators  of  the  specific  resultant  that  will 
need  to  be  employed. 

We  remark  that  the  coefficients  of  the  fi  may  involve  other  unknown  parameters.  For 
example  in  our  discussion  of  the  fire  allocation  method,  we  constructed  a  simple  example 
(2.a.ii)  that  led  to  a  system  of  four  homogeneous  quadratic  (degree  2)  equations  in  four 
variables  gi,g2,  hi,h2  whose  coefficients  involved  a  parameter  A. 

The  typical  output  of  a  resultant  depends  on  the  number  of  equations  versus  the 
number  of  variables.  It  is  important  to  note  that  when  all  the  equations  are  homogeneous, 
we  count  the  number  of  variables  as  being  one  less  than  the  actual  number.  In  this  case 
we  are  interested  in  solutions  other  than  the  trivial  solution  (all  zeros).  There  are  two 
essentially  different  cases: 

Case  1:  m  >  n  (overdetermined). 

This  is  the  case  where  we  have  more  equations  than  unknowns,  and  where  we  generally 
expect  to  have  no  solutions.  The  resultant  will  be  a  system  of  equations  (one  equation 
when  m  =  n  +  1)  in  the  symbolic  coefficients  of  the  fi,  that  has  the  following  property: 
when  we  substitute  the  specific  coefficients  of  the  /,-,  we  will  get  zero  in  every  equation  in 
the  resultant  system  if  and  only  if  the  original  overdetermined  system  has  a  solution. 

Case  2:  m  <  n  (exact  and  undetermined). 

In  this  case  the  number  of  equations  is  less  than  or  equal  to  the  number  of  variables, 
and  we  expect  to  have  solutions.  In  fact,  if  we  allow  complex  solutions  and  solutions  at 
infinity,  we  are  guaranteed  to  have  solutions. 

Of  course,  only  when  m  =  n  do  we  expect  a  finite  number  s  of  solutions.  Bezout’s 
Theorem  then  provides  a  count  oi  s  =  did2  •  •  -  dm  solutions  (counting  complex  solutions, 
solutions  at  infinity,  and  counting  with  appropriate  multiplicities).  Unfortunately,  the 
possibility  also  exists  (even  when  m  =  n),  that  there  will  be  an  infinite  number  of  solutions. 

In  general,  for  m  <  n,  the  resultant  will  be  one  equation  in  n  —  m  + 1  of  the  variables. 
In  effect,  the  resultant  makes  it  possible  for  us  to  eliminate  m  —  1  of  the  variables.  For  ex¬ 
ample,  if  we  choose  to  eliminate  Xn-m+2i  •  •  •  then  the  resultant  R  will  be  a  polynomial 
R(xi , . . . ,  Xn-m-Hi  )•  If  (ai , . . . ,  )  is  a  solution  to  i?  =  0,  then  there  will  exist  values 
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an-m+2j  •  •  • )  <^n  such  that  («!, . . . , ctn-m+i , Q:n-m+2,  •  •  •  >«n)  is  a  solution  to  the  original 
system.  (One  must  be  a  little  careful  here.  The  system  should  be  modified  to  make  it 
homogeneous  with  respect  to  Xn-m+2,-  •  • ,  by  adding  appropriate  powers  of  a  variable 
w.  The  values  a„_  to+2,  •  •  • ,  On  should  then  be  regarded  as  points  (a„_m+2  :  . .  •  :  On  :  1) 
in  projective  tn  —  1  space  .  We  must  allow  for  the  possibility  that  this  point  will  be 

at  infinity,  where  lo  =  0.  In  that  case,  a  solution  to  i?  =  0  would  not  necessarily  give  rise 
to  a  solution  of  the  original  system. 

Example:  Consider  the  system  of  m  =  2  equations  in  n  =  3  variables:  ixyz  -1  =  0  and 
y  +  xz  —  1  =  0.  The  resultant  eliminating  2  is  R{x,y)  =  x(4y^  — 4y  +  l).  When  x  =  0  we  will 
have  R  =  0,  but  clearly  our  system  has  no  solution  when  x  =  0.  However,  homogenizing 
with  respect  to  2  gives  the  system 

Axyz  —  t/;  =  0  and  (y  —  l)t4;  +  2:2  =  0. 

Now  when  we  look  at  the  condition  x  =  0  we  find  that  (2,11;)  =  (1  :  0)  is  a  solution.  This 
is  a  point  at  infinity. 

Notice  that  we  also  have  solutions  to  R  =  0  when  x  ^  0  by  taking  y  -  This 
yields  z  =  ■^.  Geometrically  the  solution  set  is  a  hyperbola  in  the  plane  y  =  1/2  in  space. 
The  resultant  “projects”  that  hyperbola  to  the  line  y  =  1/2  in  the  xy-plane,  except  that 
(x,  y)  =  (0, 1/2)  is  not  hit. 
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□ 

In  this  context  (the  underdetermined  case)  the  resultant  can  be  viewed  as  a  projection 
of  the  nominally  n  —  m  dimensional  locus  of  solutions  in  R”  to  an  n  —  m  dimensional  locus 
(hypersurface)  in  Note  that  in  our  example  n  =  3,  n  =  2,  and  we  are  projecting 

the  one- dimensional  locus  of  solutions  in  to  a  one-dimensional  locus  in  R^  which  is 
described  by  one  equation  y  —  1/2  =  0. 

Our  discussion  in  this  section  begins  with  the  first  major  distinction  in  methods, 
namely  the  one  based  on  the  number  of  variables  n.  The  case  n  =  1  of  a  single  variable  is 
discussed  first).  We  then  move  onto  the  multivariate  case  n  >  2. 

i.  Resultants  of  Polynomials  in  One  Variable 

Given  two  positive  integers  r,  s  >  1  and  two  polynomials  in  one  variable 

f{x)  =  arX^  -!-•••  +  ciix  -t-  oo  g{x)  =  b^x^  +  •  •  •  +  bix  +  bo 

of  degree  less  than  or  equal  to  r  and  s  respectively,  we  define  their  resultant  Rr,sif-i9)  by 
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Sylvester’s  formula: 

oo  <^1 

0  ao  oi 

0  0  .. 

bo  bi 
0  bo  bi 

0  0  .. 

which  is  the  determinant  of  an  r  +  s  by  r  +  s  matrix  with  s  rows  involving  the  a  s  and  r 
rows  involving  the  6’s.  (Note  that  the  a’s  and  6’s  can  be  functions  of  other  parameters, 
e.g.  A  in  our  fire  allocation  examples;  see  page  27  for  a  case  where  r  =  4,  s  =  4.) 

For  example  the  resultant  of  two  quadratic  equations  is: 

R2,2i<i2^^  +  aix  +  ao,  b2X^  +  bix  +  to)  = 

"I”  ^o^2b\  ~  aoaibib2  a^bob2  0102^0^1 
+02^0  —  20002^0^2- 

Note  that  each  monomicil  in  the  resultant  hcis  total  degree  r +s  and  that  it  is  bihomogeneous 
of  bidegree  (s,  r)  in  the  a’s  and  t’s  respectively. 

Basic  properties  of  the  resultant  Rr,sifig)’ 

1)  Relationship  to  Common  Roots 

Rr,s{fi 9)  =  o-rb's  n(^i  “  Vj)  where  xj , . . . ,  are  the  roots  of  /  and  , . . . ,  are  the 
ij 

roots  of  g.  (Here  we  are  assuming  ar  ^  0  and  bg  /  0.)  Thus  Rr,s{f,9)  will  be  zero 
if  and  only  if  /  and  g  have  a  root  in  common.  When  the  a’s  and  b’s  involve  other 
parameters,  Rr,s{f,g)  will  be  a  function  of  those  parameters.  Setting  Rr,a{f,g)  =  0 
will  then  tell  us  which  values  of  the  parameters  cause  /  and  g  to  have  a  root  in 
common. 

V 

2)  Irreducibility 

Rr,s{f-,g)  e  I[ao,...,ar,bo,...,bs]  is  irreducible,  i.e.  the  resultant  is  an  irreducible 
polynomial  with  integer  (Z)  coefficients  in  (r  +  1)('S  +  l)  =  rs  +  'S  +  r  +  l  variables. 
(Here  we  are  treating  the  coefficients  of  /  and  g  as  symbolically,  i.e.  as  independent 
variables.) 
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3)  Symmetry 

RrAf^g)  =  (--^rR.AgJ) 

4)  Factorization 

Rri+r2ysiflf2’)  9)  ^  Rri,sifl'>9)^r2,3{f2'i9)* 

ii.  Resultant  Methods  for  Systems  of  Polynomial  Equations  in  Several  Vari¬ 
ables 

The  linear  algebra  techniques  discussed  in  this  section  can  be  used  to  solve  systems 
of  polynomial  equations  in  several  variables.  If  there  are  only  two  equations,  then  the 
Sylvester  technique  (discussed  above)  can  be  employed,  by  treating  all  but  one  variable  as 
part  of  the  coefficients.  However,  when  the  number  of  equations  exceeds  two,  the  Sylvester 
approach  can  be  misleading.  For  example,  taking  the  equations  two  at  a  time  using  the 
Sylvester  determinant  can  lead  the  user  to  the  conclusion  that  there  is  a  common  solution, 
when  in  fact,  there  are  no  common  solutions  for  the  system  of  equations  taken  as  a  whole. 

What  it  means  to  “solve”  a  given  set  of  polynomial  equations  depends  upon  the  num¬ 
ber  of  variables  and  the  number  of  equations.  Assuming  the  equations  are  inhomogeneous, 
let  “n”  be  the  number  of  variables  and  “m”  be  the  number  of  equations.  The  expected 
dimensionality  of  the  set  of  solutions  is  n  m  when  viewed  over  the  complex  numbers.  For 
example,  if  there  are  three  equations  (m  =  3)  and  five  variables  (n  =  5),  then  the  space  of 
solutions  is  expected  to  have  dimension  n  —  m  =  5  —  3  =  2.  Geometrically,  the  solutions 
are  surfaces.  Sometimes,  however,  components  of  excess  dimension  occur  in  the  space 
of  solutions.  These  are  geometric  loci  of  higher  dimension  than  the  expected  dimension. 
They  occur  because,  in  a  very  loose  sense,  the  equations  have  certain  dependencies. 

Finally  a  word  about  homogeneous  equations.  Recall  that  a  set  of  polynomial  equa¬ 
tions  is  considered  homogeneous  if  in  each  equation,  all  the  terms  have  the  same  degree. 
If  this  is  not  the  case,  even  for  only  one  of  the  equations,  the  set  is  regarded  as  inhomoge¬ 
neous.  For  systems  of  homogeneous  equations,  the  number  n  of  variables  should  be  taken 
as  one  less  than  the  actual  number  of  variables  when  computing  expected  dimensions. 
The  Macaulay  Resultant  and  the  t/-Resultant 

The  Macaulay  resultant  is  the  ratio  of  two  determinants  formed  from  the  coefficients 
of  the  given  polynomials  in  a  manner  to  be  described  later  in  this  section.  If  the  number 
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of  equations  exceeds  the  number  of  variables  by  one  {n  -  m  =  -1),  then  the  Macaulay 
resultant  tests  whether  or  not  a  common  solution  exists.  (For  systems  of  homogeneous 
equations  where  the  number  of  equations  equals  the  number  of  variables,  the  expected 
dimension  is  still  -1,  and  the  Macaulay  resultant  tests  for  a  non-trivial  common  solution, 
i.e.  a  solution  other  than  (0, . . . ,  0).) 

If  there  are  as  many  inhomogeneous  equations  as  unknowns  (n  —  m  =  0),  then  the 
equations  can  often  be  solved  by  adding  the  [/-equation  (explained  later  in  this  section) 
to  the  homogenized  set  and  forming  the  Macaulay  resultant.  The  Macaulay  resultant  is 
then  called  the  [/-resultant. 

In  some  cases,  however,  there  will  be  an  excess  component  which  masks  some  or  all  of 
the  desired  solutions.  In  this  case  Canny’s  Generalized  Characteristic  Polynomial  (GCP ) 
approach  is  useful  (see  [11]). 

In  order  to  illustrate  the  various  methods,  the  following  system  of  three  polynomial 
equations  will  be  used: 

/i=y-3x -1-5  =  0 

/2  =  -|-  —  5  =  0 

fz  =  y  —  +  3x^  —  3x  -t- 1  =  0. 

Here  we  have  three  inhomogeneous  equations  in  two  variables  {n  —  m  —  2  —  3  =  — !)■ 
The  multiresultant  techniques  described  below  can  be  used  to  test  for  the  existence  of  a 
solution. 

Step  1:  Homogenization 

The  equations  must  first  be  homogenized.  This  is  done  by  adding  a  third  variable,  z. 
Specifically  x  is  replaced  by  x/z  and  y  is  replaced  by  y/z,  and  the  factors  of  z  are  cleared 
from  the  denominators.  In  the  above  example  this  leads  to  three  equations 

/i  =  y  —  3x  -|-  5z  =  0 

/2  =  x^  -I-  y^  -  5z^  =  0 

h  —  —  x^  +  3x^z  —  3xz^  -h  z^  =  0. 

This  is  the  homogenized  version  of  the  original  system. 
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step  2:  Degree  Determination 

Each  of  the  multiresultants  being  considered  involves  the  coefficients  of  various  monono- 
mials  that  appear  in  the  equations  of  the  system.  The  variables  involved  in  the  monomials 
are  the  set  of  variables  that  appear  in  the  homogeneous  form  of  the  polynomial  equations. 
For  example,  the  homogeneous  polynomial  equations  above  have  the  variables  x,  y,  and  z. 
All  the  monomials  in  a  given  equation  are  constrained  to  have  the  same  degree  because 
we  have  homogenized.  The  ^^overall  degree”  of  the  system  is  determined  from  the  degrees 
of  the  individual  homogeneous  equations  by  the  following  rule: 

m 

C?  =  1  +  ^  \{d’i  —  1) 

i=l 

where 

m  =  the  number  of  equations 

di  =  the  degree  of  the  “zth”  equation. 

For  the  homogeneous  polynomials  above  {fi,  f 2  ■,  and  )  the  degrees  are: 

Equation  Degree 

/i 

f2  d2  =  2 

/a  dz  -  3. 

Therefore, 

(i  =  l  +  (l-l)  +  (2-l)  +  (3-l)  =  4. 

(As  another  example,  in  our  simple  fire  allocation  model  we  had  four  homogeneous  quadratic 
equations  in  the  weapon  values  gi,g2,  ~  this  yields  d  =  5.) 

Step  3:  Matrix  Size  Determination 

Each  of  the  multiresultants  to  be  discussed  involves  the  ratio  of  two  determinants. 
The  numerator  is  the  determinant  of  a  matrix,  the  formation  of  which  will  be  discussed 
in  subsequent  sections.  The  denominator  determinant  is  formed  from  a  submatrix  of  the 
numerator  matrix. 

The  number  of  variables  in  the  inhomogeneous  equations  is  n.  Since  one  additional 
variable  has  to  be  added  to  homogenize  the  equations,  the  number  of  variables  in  the 
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homogeneous  equations  is  n  +  1.  The  size  of  the  numerator  matrix  equals  the  number 
of  monomials  in  the  n  +  1  variables  that  have  overall  degree  d  (discussed  in  the  previous 
section). 

Numerator  Matrix  Size  = 

For  the  three  polynomial  equations  (/i,/2,/3)  we  have  already  calculated  that  d  =  4. 
Since  the  original  set  of  inhomogeneous  variables  consisted  of  x  and  y,  n  equals  2.  Thus 
for  our  example, 

Numerator  Matrix  Size  = 

i.e.  it  is  a  15  X  15  matrix.  (In  the  fire  allocation  case  mentioned  above  the  numerator  size 
would  be  (5)  =  56  which  is  rather  large.) 

Step  4:  Determining  “Big”  vs.  “Small”  Exponents 

A  few  of  the  15  monomials  involving  the  variables  x,y,  and  z  with  an  overall  degree 
of  4,  include: 

(1) 

(2)  ar^y^. 

In  the  next  step,  we  will  discuss  whether  certain  monomials  are  reduced.  This  will  be 
determined  by  whether  the  exponents  are  “big”  or  “small”.  In  this  step  we  discuss  how 
“bigness”  is  defined. 

Each  variable  will  be  associated  with  a  particular  equation.  For  example  the  first 
variable,  x,  will  be  associated  with  the  first  equation,  /i.  The  second  variable,  y,  will 
be  associated  with  the  second  equation,  /2,  etc.  The  degrees  of  the  associated  equations 
define  “bigness”  for  the  exponents  of  that  variable.  Specifically,  since  di  (the  degree  of  fi ) 
is  1,  if  the  exponent  of  x  is  greater  than  or  equal  to  1,  it  is  considered  big.  Since  ^2  =  2, 
whenever  the  exponent  of  y  is  greater  than  or  equal  to  2,  it  is  considered  big.  The  degree  of 
/a  is  3,  therefore,  whenever  the  exponent  of  z  is  greater  than  or  equal  to  3,  it  is  considered 
big. 

For  example,  consider  the  monomial  numbered  (1),  namely,  yz^.  The  exponent  of  y 
is  1.  This  is  less  than  da,  and  is  considered  small.  The  exponent  of  z  is  3.  This  is  equal  to 


2  +  4' 
4 


6! 


(2!)(4!) 


=  15 
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dz,  and  is  therefore  big.  On  the  other  hand,  if  we  consider  the  monomial  numbered  (2), 
namely,  x^y^.  The  exponent  of  x  is  2.  This  is  greater  than  dj,  and  is  big.  The  exponent 
of  y  is  2.  This  is  equal  to  <^2,  and  is  big. 

Step  5:  Determining  the  Reduced  Monomials 

If  for  a  particular  monomial  of  degree  d  the  exponent  of  only  one  variable  is  big,  the 
monomial  is  said  to  be  reduced.  In  the  previous  step  only  monomial  (1)  is  reduced.  For 
monomial  (1)  only  the  exponent  of  z  is  big;  whereas  for  (2),  both  the  exponent  of  x  and 
the  exponent  of  y  are  big.  Thus  monomial  (2)  is  not  reduced. 

The  Macaulay  Resultant 

The  Macaulay  Resultant  is  the  ratio  of  two  determinants.  The  numerator  is  the  deter¬ 
minant  of  a  matrix  which  we  will  call  the  A  matrix.  The  denominator  is  the  determinant 
of  a  matrix  which  we  will  call  the  M  matrix 

det  \A\ 
detlMf 

Step  6:  Creating  the  A  Matrix: 

We  have  discussed  above  how  the  size  of  the  A  matrix  is  determined.  In  this  step  we 
will  show  how  the  matrix  elements  are  obtained. 

Each  row  and  column  of  the  matrix  should  be  thought  of  as  being  labeled  by  one  of 
the  monomials  of  degree  d.  Recall  that  for  /i,/2,  and  /s  there  were  15  possible  monomials 
of  degree  4  in  x,  y,  z,  and  therefore  the  A  matrix  would  be  15  x  15. 

There  are  three  rules  for  determining  the  elements  of  the  A  matrix.  After  presenting 
the  rules,  the  example  involving  /i,/2,  a.nd  /j,  will  be  used  to  illustrate  the  process.  The 
reader  may  find  it  helpful  to  read  the  example  simultaneously  with  the  rules. 

Rules  for  the  elements  of  each  column; 

(1)  Search  the  monomial  lal^eling  that  column  from  left  to  right  for  the  first  variable 
with  a  big  exponent.  Such  a  variable  must  exist.  Call  it  the  marker  variable. 

(2)  Form  a  new  polynomial  from  the  polynomial  associated  with  this  marker  variable 
by  multiplying  the  associated  polynomial  by  the  monomial  and  dividing  by  the  marker 
variable  raised  to  the  degree  of  the  associated  polynomial. 
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(3)  The  coefficients  of  the  new  polynomial  are  the  elements  of  the  columns.  Each 
coefficients  goes  in  the  row  labeled  by  the  monomial  it  multiplies.  All  the  other  rows  get 
zeroes. 

Recall  that  for  our  example  system  of  equations  /i,/2,  /a  there  are  15  monomials  of 
degree  4  that  can  be  formed  from  x,  y,  and  2.  Two  of  these  were  considered  above,  namely 
yz^  and  x^y^. 

(a)  For  the  column  labeled  by  yz^: 

(1)  The  first  variable  with  a  big  exponent  is  2,  so  2  is  the  marker  variable. 

(2)  The  polynomial  associated  with  2  is  /a.  Multiply  /a  by  the  monomial  y2^  and 
divide  this  product  by  2^. 


-  x^y  +  Zx‘‘yz  -  Zxyz^  +  yz\ 


(3)  The  coefficient  of  y'^z'^  is  +1.  Therefore  the  element  of  the  row  labeled  y^z  is 
+1.  The  coefficient  of  x^y  is  -1.  Therefore  the  element  of  the  row  labeled  x^y  is 


-1.  The  coefficient  of  x'^yz  is  +3.  Therefore  the  element  of  the  row  labeled  x^yz 
is  +3.  The  coefficient  of  xyz'^  is  -3.  Therefore  the  element  of  the  row  labeled 
xyz^  is  —3.  The  coefficient  of  yz^  is  +1.  Therefore  the  element  of  the  row  labeled 
yz^  is  +1.  All  other  entries  in  the  column  are  zero. 

(b)  For  the  column  labeled  by  x^y^: 

(1)  The  first  variable  with  a  big  exponent  is  x,  so  x  is  the  marker  variable. 

(2)  The  polynomial  associated  with  x  is  /i.  Multiply  /i  by  the  monomial  x^y^,  and 
divide  this  product  by  x. 


X  X 

(3)  The  coefficient  of  xy^  is  +1.  Therefore  the  element  of  the  row  labeled  xy®  is  +1. 
The  coefficient  of  x^y^  is  -3.  Therefore  the  element  of  the  row  labeled  x^y^  is 
-3.  The  coefficient  of  xy^2  is  +5.  Therefore  the  element  of  the  row  labeled  xy^2 
is  +5. 
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When  all  the  columns  are  determined,  the  A  matrix  in  our  example  takes  the  form: 

A  MATRIX 


x^ 

x^ 

x^ 

x2 

x2 

X^ 

X 

X 

X 

X 

y 

y^ 

y 

y’ 

y^ 

y 

y^ 

y^ 

y^ 

y 

2 

2 

22 

2 

22 

2^ 

2 

..3 

2^ 

-3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

x^ 

y 

1 

-3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

x3 

2 

5 

0 

-3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

x2 

y" 

0 

1 

0 

-3 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

X2 

y 

2 

0 

5 

1 

0 

-3 

0 

0 

0 

0 

0 

'  0 

1 

0 

3 

0 

22 

0 

0 

5 

0 

0 

-3 

0 

0 

0 

0 

0 

0 

1 

0 

3 

X 

y" 

0 

0 

0 

1 

0 

0 

-3 

0 

0 

0 

0 

0 

0 

0 

0 

X 

y^ 

2 

0 

0 

0 

5 

1 

0 

0 

-3 

0 

0 

0 

0 

0 

0 

0 

X 

y 

22 

0 

0 

0 

0 

5 

1 

0 

0 

-3 

0 

0 

0 

0 

-3 

0 

X 

2^ 

0 

0 

0 

0 

0 

5 

0 

0 

0 

-3 

0 

0 

0 

0 

-3 

y" 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

y" 

2 

0 

0 

0 

0 

0 

0 

5 

1 

0 

0 

0 

1 

0 

0 

0 

y^ 

22 

0 

0 

0 

0 

0 

0 

0 

5 

1 

0 

-5 

0 

1 

1 

0 

y 

2^ 

0 

0 

0 

0 

0 

0 

0 

0 

5 

1 

0 

-5 

0 

1 

1 

2^^ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

5 

0 

0 

-5 

0 

1 

The 

determinant  of 

the  above  A  matrix 

is  zero. 

If  the  M  matrix  is  nonzero,  this  would 

imply  that  the  system  has  a  solution. 


Step  7;  Creating  the  M  Matrix: 

The  denominator  of  the  Macaulay  Resultant  is  the  determinant  of  the  M  matrix.  The 
M  matrix  is  a  submatrix  of  the  A  matrix.  It  consists  of  the  elements  which  have  row  and 
column  monomial  labels  which  are  not  reduced.  Recall  that  a  monomial  is  not  reduced  if 
it  has  more  than  one  variable  with  a  big  exponent. 

The  size  of  the  M  matrix  equals  the  size  of  the  A  matrix  minus  D,  where 

^  m 

i=l  ijij 

In  our  example, 

D  =  d2d,  +  dids  +  did2  =  (2)(3)  +  (1)(3)  +  (1)(2)  =  11. 


40 


So  that  the  size  of  the  M  matrix  is  15  -  11  =  4.  The  actual  M  matrix  for  /i,  /2  and  /s  is: 


M 

MATRIX 

2  2 
^  y 

xy* 

xy'^z 

X2^ 

x^y^ 

-3 

0 

0 

0 

xy^ 

1 

-3 

0 

0 

xy^z 

5 

0 

-3 

0 

XZ 

0 

0 

0 

-3 

The  determinant  of  this  M  matrix  yields  a  value  of  81.  Since  the  A  matrix  was  zero, 
the  Macaulay  Resultant  is  zero  which  implies  that  there  is  a  solution  to  our  system.  A 
quick  plot  of  the  three  polynomials  (/i ,  /2,  and  /a)  confirms  that  there  is  a  common  point 
at  x  =  2  and  y  =  1. 

Sometimes  both  the  A  matrix  and  the  M  matrix  have  zero  determinant.  This  inde¬ 
terminacy  can  often  be  circumvented  if  the  polynomials  are  first  written  with  symbolic 
coefficients.  The  determinants  of  the  A  and  M  matrices  are  then  computed,  and  polyno¬ 
mial  division  is  performed.  At  that  point  the  symbolic  coefficients  are  replaced  by  their 
numerical  values  to  check  if  the  resultant  is  zero.  Since  one  does  not  know  ahead  of 
time  whether  or  not  this  “division  by  zero”  condition  will  arise,  the  symbolic  coefficient 
approach  is  the  best  strategy.  It  is  also  often  sufficient  to  treat  just  a  subset  of  the  co¬ 
efficients  symbolically  -  sometimes  as  few  as  a  single  symbolic  coefficient  will  remove  the 
indeterminacy. 

The  17-Resultant 

For  problems  with  as  many  inhomogeneous  equations  as  variables,  the  U  resultant  can 
often  be  used  to  solve  for  the  point  solutions.  The  three  polynomial  equations  /i ,  /2,  /a  do 
not  satisfy  these  conditions  since  there  are  three  equations  in  two  inhomogeneous  variables, 
X  and  y.  However,  if  we  take  just  the  first  two  equations,  namely  /i  and  /2,  we  would  have 
a  system  with  as  many  equations  as  variables. 

The  given  equations  must  first  be  homogenized.  This  adds  one  additional  variable.  We 
then  add  one  additional  equation  to  the  system.  This  equation  is  called  the  [/-equation. 
If  X  and  y  are  the  given  variables  and  2  is  the  homogenizing  variable,  then  the  U  equation 
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tcikes  the  form: 

uix  +  U2y  +  U3Z  =  0. 

The  Macaulay  Resultant,  i?,  is  then  computed  for  these  m  + 1  equations,  treating  the 
Ui  as  symbolic  coefficients.  The  result  is  called  the  [/-resultant.  Notice  that  R  will  be  a 
polynomial  in  the  Ui's  and  the  coefficients  of  the  original  equations. 

After  R  is  determined,  it  is  then  factored  into  linear  factors.  For  each  linear  factor 
there  is  a  point  solution  of  the  original  system  of  equations.  The  coordinates  of  each 
solution  are  determined  as  ratios  of  the  coefficients  of  the  u’s.  The  denominator  is  always 
the  coefficient  of  the  Ui  associated  with  the  homogenizing  variable.  In  our  example,  this  is 
the  coefficient  of  U3.  Thus 

coeff.  of  ui  _  coeff.  ofu2 
coeff.  of  U3  ^  coeff.  of  U3 

For  example,  if  a  linear  factor  turned  out  to  be 


Ui  —U2  —  «3, 


then  the  coordinates  of  the  associated  solution  would  be 


As  mentioned  above,  we  can  use  the  [/-resultant  to  solve  fi  and  /2  for  x  and  y  {n—m  = 
0).  In  this  example,  we  will  also  demonstrate  the  symbolic  approach  alluded  to  in  the 
previous  section.  Recall  that  the  homogenized  form  of  f\  and  /2  is: 

=  y  —  Sr  -b  52  =  0 
f2=x^+y^-^z^  =  Q. 

Rewriting  these  two  equations  with  symbolic  coefficients  and  including  the  [/-equation 
yields: 

/i  =  Oil  -I-  hiy  -b  Ci2  =  0 
/z  —  0-2^^  "b  "b  ^22^  =  0 
U  =  till  -b  U2y  +  U32  =  0 
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where,  ai  =  —3;  bi  =  1;  Cj  =  5;  02  —  1;  ^2  —  1;  ^2  —  —5. 

The  [/-resultant  is  calculated  in  the  same  way  as  the  Macaulay  resultant,  i.e.  with  the 
A  matrix  and  the  M  matrix,  except  now  we  are  using  symbolic  coefficients. 


A  MATRIX 

xy 

xz 

y^ 

yz 

^2 

Oi 

0 

0 

0-2 

0 

0 

xy 

h 

Oi 

0 

0 

Ui 

0 

xz 

Cl 

0 

ai 

0 

0 

Ui 

y^ 

0 

h 

0 

62 

U2 

0 

yz 

0 

Cl 

hi 

0 

U3 

U2 

0 

0 

Cl 

C2 

0 

U3 

The  corresponding  M  matrix 

is  a 

single  element,  namely  oi. 

The  determinant  of  M  is  divided  into  the  determinant  of  A  to  obtain  the  [/-resultant. 
Finally,  the  symbolic  coefficients  are  replaced  by  their  numeric  equivalents.  (This  could 
have  been  done  from  the  outset,  unless  gi  had  been  zero.)  The  result  is 


10(ui  —  2u2  +  U3)(2ui  -1-^2+  Us)- 

This  yields  two  solutions: 

Solution  ^1: 


coeff.  of  ui  -t-1  ,  ,  coeff.  of  U2 

- —  =  —  =  -f-1  and  y  =  - - — z — 

coeff.  of  U3  -1  coeff.  of  U3 


-2 

+I 


=  -2. 


Solution  ^2: 


coeff.  of  ui  +2  ^  ,  coeff.  of  U2 

- =  —  =  +2  and  v  = - := — :: — 

coeff.  of  U3  -|-1  coeff.  of  U3 


We  remark  that  the  [/-resultant  will  be  identically  zero  and  give  no  information  if  the  set 
of  common  solutions  contains  a  component  of  excess  dimension  one  or  more.  Moreover, 
this  component  may  be  at  infinity  where  the  homogenizing  variable  is  zero. 

The  GCP  Approach 

The  Generalized  Characteristic  Polynomial  (GCP)  approach  avoids  the  problem  of 
components  of  excess  dimension  in  the  space  of  solutions.  It  can  be  used  together  with 
the  [/-resultant  which  was  discussed  above.  If  the  [/-resultant  leads  to  an  indeterminant 
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(0/0)  form  even  when  symbolic  coefficients  are  used,  an  ^‘excess”  solution  exists.  The  GCP 
takes  the  form 

R  =  evaluated  at  s  =  0  (after  division) 

det  \M  -  sl\ 

where  A  and  M  are  the  matrices  defined  above,  s  is  a  perturbation  parameter,  and  I  is 
the  identity  matrix. 

One  way  to  carry  out  the  above  operation  is  the  following; 

(1)  Set  up  the  A  matrix  (as  described  above).  Subtract  an  s  along  the  diagonal. 
Evaluate  the  determinant.  Retain  the  coefficient  of  the  lowest  surviving  power  of 
s. 

(2)  Repeat  (1)  for  the  M  matrix. 

(3)  Divide  the  result  of  (1)  by  the  result  of  (2). 

All  of  these  multiresultant  techniques  have  one  thing  in  common.  They  require  that 
there  be  one  more  equation  than  variable,  n  —  m  =  —  1.  If  there  are  as  many  equations 
as  variables  n  -  m  =  0,  then  the  U  equation  is  added  and  the  effective  situation  is  again 
n  —  m  =  —  1.  If  there  are  more  variables  than  equations  (n  —  m  is  a  positive  integer), 
then  enough  of  these  variables  must  be  regarded  as  parameters  in  the  coefficients,  so  that 
effectively  n  —  m  =  —  1.  Geometrically  this  amounts  to  projecting  the  locus  of  solutions  to 
a  hypersurface  in  a  lower  dimensional  Euclidean  space.  Finally,  if  the  number  of  equations 
exceeds  the  number  of  variables  by  more  than  one  {n  —  m  <  —2),  then  some  technique 
other  than  the  above  Macaulay  resultant  (e.g.  a  systems  of  resultants)  must  be  employed 
to  determine  if  a  solution  exists. 

c.  Similarities  with  and  Differences  from  other  Methods. 

The  fire  allocation  method  is  by  far  the  most  general  method  we  discuss  in  this  paper. 
It  couples  fire  allocation,  attrition  or  kill  rates,  and  weapon  values  in  a  natural  way,  and 
leads  to  a  system  of  non-linear  equations  in  the  weapon  values  and  the  proportionality 
constant  A  which  has  units  of  (time)“^. 

This  system  can  in  turn  be  written  as  two  systems  of  non-linear  equations  -  one  for 
the  blue  values  and  one  for  the  red  values  -  linked  by  the  common  “eigenvalue”  A^.  The 
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similarity  of  this  result  with  the  simultaneous  eigenvalue  problem  arising  in  the  eigenvalue 
or  potential/antipotential  method  is  quite  striking.  Unfortunately,  we  have  not  yet  been 
able  to  prove  an  existence  or  uniqueness  result  for  these  non-linear  systems.  In  the  poten¬ 
tial/antipotential  case,  the  existence  and  uniqueness  is  assured  under  some  mild  hypothe¬ 
ses.  We  have  however  done  some  numerical  experiments  that  have  led  us  to  conjecture 
that  these  non-linear  systems  do  indeed  have  solutions  and  that  they  are  unique. 

The  major  drawback  with  the  fire  allocation  models  discussed  above  is  that  solving 
coupled  system  of  non-linear  equations  is  not  an  easy  task.  The  resultant  methods  we 
suggest  grow  exponentially  in  size  with  the  number  of  weapon  types.  Dealing  with  50  plus 
types  per  side  would  require  very  large  computations.  Perhaps  some  blend  of  symbolic 
and  numeric  techniques  could  be  employed  to  make  such  a  computation.  Canny  [11]  has 
successfully  explored  such  an  approach  in  computing  very  large  resultants. 

Another  drawback  is  that  the  pure  kill  rates  (the  dji  and  Cij  that  appear  in  our 
example)  need  to  be  known  in  advance.  The  best  approach  here  is  to  make  use  of  the 
results  of  a  high  resolution  low  level  simulation  (e.g.  COSAGE)  where  fire  allocation  can 
be  more  or  less  determined,  so  that  along  with  the  killer-victim  scoreboards,  the  pure 
attrition  coefficients  can  be  estimated. 

Finally,  on  the  positive  side,  we  believe  that  the  iterative  scheme  we  have  outlined 
for  coupling  attrition,  fire  allocation,  and  weapon  values  is  potentially  quite  useful.  It 
certainly  warrants  further  investigation. 
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3.  Weapon  Importances  in  ATCAL 
a.  The  Methods 

ATCAL  is  a  calibrated  attrition  model  used  to  determine  attrition  in  a  theater  level 
simulation  called  CEM  at  the  U.S.  Army  Concepts  Analysis  Agency.  Certain  parameters 
are  determined  (calibrated)  using  the  results  of  a  high  resolution,  low  level  (division  level 
for  example)  simulation  called  COSAGE.  Multiple  COSAGE  runs  are  done  with  a  mix  of 
all  the  weapons  systems  available  in  order  to  capture  all  possible  interactions. 

ATCAL  weapons  importances  are  calculated  from  entries  in  the  killer-victim  score- 


board.  Let 


(AMj)i 

Mj/^T 


3>lici  — 


(AiV,),- 

NiAT 


for  1  <  z  <  n  and  1  <  j  <  m.  Here  {AMj)i  is  the  number  of  type  j  killed  by  opposing 
type  i  during  the  engagement,  AT  is  the  length  of  the  engagement,  and  Mj  is  the  initial 
number  of  type  j  present  at  the  beginning  of  the  engagement.  Thus  aij  is  the  fraction  of 
type  j  killed  by  opposing  type  i  per  unit  time.  Likewise  bji  is  the  fraction  of  type  i  killed 

n 

by  opposing  type  j  per  unit  time.  Note  that  ^  akj  equals  the  fraction  of  j  killed  per  unit 

k-i 


m 

time  and  ^  hsi  equals  the  fraction  of  i  killed  per  unit  time. 

3=1 


The  weapon  importances  Aj ,  z  =  1 , . . . ,  n  and  Yj,  j  —  1 , . . . ,  m  are  defined  by: 


m 

E 

n  ^3 

for  z  =  1, . . . ,  n 

jfc=l 

n 

E 

{haf  ^ 

for  j  =  1,  ...  ,771 

i=l 

5=1 

Xi  >  0 


Yj>0. 


Notice  that  =  0  for  all  z  and  Yj  =  0  for  all  j  is  a  solution  to  the  above  system  of 
equations.  Clearly  we  want  a  non-trivial  solution  where  not  all  these  values  are  zero.  Note 
also  the  Xi  is  the  total  importance  of  all  weapons  of  that  type  as  is  Yj. 
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In  this  system  of  non-linear  equations  the  quantity 


m 


E*'.. 


AJVi 


is  the  fraction  of  all  kills  of  i  that  were  caused  by  j.  Likewise 


dij 

n 

k^\ 


AMj 


is  the  fraction  of  all  kills  of  j  that  were  due  to  i. 

We  will  discuss  the  motivations  advanced  for  these  equations  in  section  3c.  For  now 
we  content  ourselves  with  showing  one  important  property  exhibited  by  this  formulation 
of  weapon  importances,  namely  the  so-called  divisibility  property.  It  says  that  if  we  divide 
a  type  into  two  types  the  total  importance  assigned  to  them  does  not  change. 

To  prove  this,  suppose  we  divide  Ni  for  type  i  =  1  into  N[  -(-  N"  =  Ni  and  apportion 


losses  and  kills  accordingly: 


iV:  ^  iVi  ^ 

,  r(AMj)i 

=  (1  -  ^)«ij 
,f  r(AiVi)i  , 

o,-i  =  r  w.  ^  —  Oji 


(riVi)AT 

u" 

hi 

Our  new  set  of  importance  equations  are 


b'U  =bji. 


j=l  ^  Qkj  j=l  Clkj 

k=l 


k=l 

m  /  n3 


(«iir 


m  (  II  \3 

j=l  ^  Ukj  >=1  13  ^kj 

k=l  k=l 

^  I'  \3 

J~l  13  ^kj 

k=l 


(1  -  r)^i3- 
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j  =  1,. . .  ,m. 


(Note:  a[j  +  a'(j  +  a2j  +  •  •  •  +  Qnj  =  aij  -\ - +  Onj-) 


y?  = 

J  m 


X'  -j- 

E  t:.  ‘  E  E  t-i 

j=:X  S — 1  <3 — 1 


(Note:  <5k>l  .  =  ^.) 

E*;.  E‘v.  E‘-' 


8=1  8=1 


Now  if  (Xi , . . . , X„,  Fi, . . . ,  Fm)  solves  the  original  system  then 

(x(,xi',X2,...,x„,yi,...,y„.)  =  (rXi,(i-r)Xi,X2,...,x„,yi,...,y^) 

solves  the  new  system.  Conversely,  if  {X[,  X'^ ,  X2,  •  •  • ,  Fm)  satisfies  the  new  system  then 
X'  X'' 


1— r 


and 


(Xi,X2,...,F,„)  =  |^^,X2,...,F,, 

satisfies  the  original  equations.  In  particular 

X'{ 


Xi 


1  —  r 


so  that 

x[  +  X';  =  rXi  +  (1  -  r)Xi  =  Xi. 

This  means  that  the  total  value  of  the  prime  and  double  prime  weapon  types  equals  the 
original  total  value  Xi  for  type  i  =  1.  This  proves  the  desired  divisibility  property. 

Notice  that  there  is  nothing  sacred  in  this  proof  about  using  cubes.  Any  power  a 
works  as  long  as  the  exponents  match 

fc=l 

=  tor  i  = 

«=1  £  ^si 

S=1 

This  makes  the  choice  of  cubes  seem  somewhat  arbitrary.  We  remark  that  the  existence 
and  uniqueness  proof  below  for  a  =  3  will  also  work  for  arbitrary  o. 
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b.  Existence  and  Uniqueness  of  Solutions  to  the  ATCAL  Weapon  Importance 
Equations. 

ATCAL  has  been  criticized  on  a  number  of  counts  (see  Anderson  [56]),  including  the 
lack  of  any  proof  that  its  weapon  importance  equations  always  have  a  solution  and  that 
such  a  solution  is  unique.  In  this  section  we  provide  that  proof.  While  this  is  largely  a 
mathematical  exercise,  the  proof  yields  two  important  insights.  The  first  is  an  efficient 
algorithm,  based  on  an  iterative  fixed  point  method,  for  computing  the  importances,  and 
the  second  is  the  observation  that  there  is  a  striking  similarity  between  this  proof  and  the 
fixed  point  proof  of  the  Perron- Frobenius  Theorem  (see  §lc)  which  is  the  key  result  in  es¬ 
tablishing  existence  and  uniqueness  for  the  so-called  eigenvalue  or  potential/ antipotential 
method  of  computing  weapon  scores.  This  observation  leads  to  speculation  that  the  eigen¬ 
value  scores  and  the  ATCAL  importances  will  suffer  from  similar  defects. 

We  must  establish  our  existence  and  uniqueness  result  for  a  system  of  polynomial 
equations  of  the  form: 

n 

(1)  Xf  ='^bijXj  i  =  l,...,n 

j=i 


where  n  >  2  and  bij  >  0  for  every  1  <  i,j  <  n.  We  seek  a  solution  x  — 
—  {(0, . . . ,  0)^},  that  is,  a  solution  with  all  Xi  >0  but  not  all  Xi  =  0.  Let 


9:  R|o  -  {(0, .  • . ,  0)^}  ^  R|o  -  {(0,  •  •  • ,  0)^} 


Xi 


be  the  mapping  which  sends  x  = 


to  6(x)  = 


1/3 


In  order  to  show  uniqueness,  we  must  make  one  assumption, 
that  the  null  space  of  the  matrix  B  =  {bij)  intersects  the  positive 
the  zero  vector  (0,  ...,0)^.  This  means  that  Bx  £  R>o  —  {(0, ... 
R"  g  —  {(0, . . . ,  0)^}  and  is  a  very  mild  assumption  given  that  all 
non-negative. 


Namely,  we  assume 
orthant  R>o  only  in 
,0)^}  for  every  x  € 
the  entries  of  B  are 
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Now  consider  the  standard  (n  —  l)-simplex  A”  ^  =  {(^i?  •  •  •  ^  ^>o  such  that 

^  Xj  =  1,  and  Xi  >  0  all  i  =  1, . . . ,  n}  in  R^g 

t=i 

A”-i  — ^  A""^ 

Xi 


X  = 


.Xr 


d{Bx) 

a{diBx)) 


where  (t:  R"  ^  R  sends  (xi, . . . ,  x„)^  to  ^  x^.  By  the  Brouwer  fixed  point  theorem,  there 

1=1 

is  a  vector 

'u;i 

=  I  •  1  ^  A"“^ 


with 


Setting  X  =  cr{9(Bwo)),  we  have 


.Wr 


Wo 


OjBwo) 
a{9{Bwo)) 

XWo  =  9{Bwo) 


for  a  value  X  >  0-  Equivalently, 


( w\\  fWx 

=  B 


Wn. 


.Wr 


This  is  almost,  but  not  quite,  a  solution  to  (1).  If  we  scale  by  A,  we  observe  that 


\{\Wnf 


I 

=  x^^^  I  :  I  =A^B 


Wi 


.  Wr 


=  yB 


Aiui 


. \Wn  . 


If  we  now  pick  A  =  then 


vo- 


\wi 


.\Wn. 


=  Ar^o 
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solves  (1),  i.e. 


This  establishes  the  existence  of  a  non- negative,  not  identically  zero  solution  to  (1). 
We  now  turn  to  uniqueness.  By  the  argument  above  we  know  that 

n 

xf  =J2bijXj  bij>0  i  =  l,...,n 
i=i 

has  at  least  one  solution  (under  a  mild  hypothesis  on  B  =  (bij)),  call  it 


Vl 


Vo  = 


We  will  assume  all  Vi  >  0.  (For  example  one  can  throw  out  all  weapons  with  zero  value 
and  the  corresponding  equations.) 

Now  set  Yi  =  ^  to  get  a  new  system  of  equations 


(2) 


j=l 


where  Cij  =  >  0  and  =  1  for  alH  =  1, . . . ,  n.  This  new  system  has  I  j  as  a 

j=i  '  ' 

solution. 

Suppose  there  is  another  solution  to  (2) 

/  yi 


Vi  >0 


\yn, 


with  some  t/i  >  1.  Let  io  be  the  index  for  which  yi^  is  a  maximum.  We  have 


fio  yio  ~  ^^^iojyj  I  ]  y«o 

3=1  \j=l 

n 

j=i 
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but  -  y,o  >  0  while  Ci^jiyj  -  yio)  <  0-  This  yields  a  contradiction.  Thus  if 

\ynj 

is  another  solution,  then  all  yj  are  less  than  or  equal  to  one  and  at  least  one  yi  is  less  than 
one.  (We  continue  to  assume  all  yi  >  0.) 

Let  io  be  the  index  for  which  y^g  is  a  minimum.  Then 

n 

Vio  “  2/to  ^  ^  Vj) 

n 

but  yi^  -  >  0  while  Y.  -yj)<  0- 

i=i 

We  caution  that  this  “proof’  of  uniqueness  is,  strictly  speaking,  incomplete.  What 
we  have  shown  is  that  we  can’t  have  two  different  solutions  with  all  weapon  importances 
being  positive.  But  we  can’t  rule  out  multiple  solutions  when  some  yi  are  zero. 

This  problem  with  uniqueness  breaking  down  is  virtually  identical  to  what  happens 
in  the  eigenvalue  method  when  multiple  solutions  occur.  Moreover,  our  existence  proof 
is  virtually  identical  to  the  existence  proof  for  the  eigenvalue  method  given  in  Section  1. 
One  is  tempted  to  conclude  that  both  the  eigenvalue  method  and  the  ATCAL  method  will 
exhibit,  at  least  qualitatively,  similar  behaviour.  In  particular  both  methods  will  suffer 
from  similar  defects  in  degenerate  cases,  and  any  criticism  of  the  eigenvalue  method  is 
likely  to  apply  equally  to  the  method  in  ATCAL. 

Finally  the  existence  proof  above  provides  a  very  stable  iterative  algorithm  for  solving 
the  importance  equations.  Specifically,  compute  the  {i  +  l)st  iterate  by 

e{Bxi) 

'  Xi-i-i  = - z: — 

a{9{Bxi)) 

starting  with  =  (n’  •  •  •  ’  n)  fixed  point  is  found,  we  scale  as  indicated  in 

the  proof  to  get  a  solution. 

c.  Problems  and  Similarity  with  the  Eigenvalue  Method. 


52 


The  ATCAL  weapon  importance  equations  have  a  somewhat  ad  hoc  feel  to  them.  We 
were  not  able  to  find  any  derivation  of  them  from  first  principles,  as  we  can  for  the  other 
methods.  Moreover,  as  indicated  in  3a,  the  cubic  exponent  is  not  the  only  exponent  that 

will  produce  a  satisfactory  notion  of  value. 

On  the  positive  side,  the  importance  equations  reduce  to  the  eigenvalue  method  when 

there  is  only  one  weapon  per  side: 

Xf  =  0^2X2 

=  aliXi 


where  ai2  =  ^  and  021  =  ^  are  the  fractional  amounts  killed 

per  unit  time.  The  solution  to  these  equations  is 

Y  _  3/4  1/4 
—  ^12  ^21 

V  _ 

^2  —  ^12  ^21  * 


Thus  the  value  of  a  single  unit  of  type  1  is  Vi  =  Xi/Ni  and  of  type  2  is  V2  =  X2IN2. 
Normalizing  Vi  to  1  gives 

V2  _  X2  ^  /Q21 

Vi  N2  Xi  N2  V  ai2 

Conversely  the  eigenvalue  method  for  two  weapons,  one  Red  and  one  Blue,  gives 


CBxVB  =  kVR 
CRxVR  =  iVR 


{AR)b  ^  XR 
B- AT  BAT 
_  (AB)fi  _  AB 
R- AT  RAT 


where  CB  =  CR  =  s/X  and  where  we  normalize  VB  to  1.  The  simultaneous  eigenvalue 
equations  are 

XV B  =  M  VB 
XV R  =  U  VR 

and  we  see  that  X  =  k£  =  aRBO.B-R  where 


AB  AR 

aRB  =  aBR  - 


We  thus  also  get 


^  _  5  laRB 
y  k  R  \  asR 
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which  is  the  same  as  ATCAL  gave  above  (with  indices  1  and  2  instead  of  B  and  R). 

One  odd  feature  of  the  ATCAL  importances  is  that  they  depend  on  the  units  of 
time  used.  The  unit  of  importance  is  (tirne)"^  which  is  an  intensity  (by  analogy  with 
Lanchester  and/or  the  eigenvalue  method).  However,  for  relative  importances  this  effect 
cancels  out.  Moreover  there  is  a  certain  logical  inconsistency,  in  that  the  ATCAL  attrition 
model  assumes  exponential  averaging,  i.e.  that  attrition  follows  an  exponentially  decaying 
curve  over  time,  while  the  importance  equation  are  based  on  constant  fractional  loss  per 
unit  time.  Why  not  incorporate  the  exponential  decay  into  the  model?  In  fact,  one  could 
compute  fractional  losses  as  a  function  of  time  using  COSAGE  data.  This  would  give  time 
dependent  importances  which  could  be  integrated  over  the  duration  of  the  engagement  to 
give  dimensionless  numbers  for  weapon  values. 

One  possible  derivation  of  the  ATCAL  importance  equations  is  given  below,  but  it  is 
only  heuristic,  and  is  based  on  very  poor  assumptions  about  times  to  completion  of  various 
interactions.  Recall  that  in  the  eigenvalue  method; 

n 

CB  xVBi  =  Y^  KijVRj 

i=i 

m 

CR  X  VRj  =  ^  LjiVBi 

i=l 

where  CB  =  CR  —  \/A  =  I  (intensity).  In  other  words,  a  weapon  s  value  is  proportional 
to  the  rate  at  which  it  destroys  opposing  value  (at  a  given  instant  in  time).  This  leads  to 

CB  X  yS,  X  Bj  =  ^  X  VR,  X  Rj 

j=\  ^ 

or 

CB  X  TVs,  =  Y.  ^ 

j-1  ^ 

where  TV’s,,  is  the  total  value  of  Rlue  weapon’s  of  type  i,  etc.  This  yields: 

and 

TV,,^±h,.TVs.[^) 

i=l  ^  ^ 
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where  the  units  on  aij  and  hji  are  (time)  ^  and  the  units  on  CB  —  CR  —  —  I  are  also 

(time)-^  Now  regard  TFb,.  =  Xs.dt  =  where  Xg,  is  the  total  value  added 

to  type  Bi  per  unit  time.  Here  we  assume  that  Xb^  is  constant  over  a  short  time  interval 
AT.  Note  Xb;  has  units  (time)“^. 

The  short  time  rate  of  contribution  to  value  in  the  engagement,  i.e.  importance,  is  then 
proportional  to  the  rate  of  contribution  to  opposing  total  value  that  the  weapon  eliminates: 


Xb.  =Y^a,jXR.  . 


Now  relax  assumption  that  CB  is  constant,  making  it  dependent  on  the  i^j  pair.  This 
is  a  kind  of  disaggregation  into  separate  interactions!  Set  '"'‘CBi^j"  =  yjXsi Xr-  which  has 
units  =  (time)“^  and  view  '^CBij'  as  a  time  constant  for  the  i,j  interaction.  If  we  then 
set  fxj  equal  to  the  fractional  involvement  of  Bj  in  kills  of  Rj,  i.e.  fij  =  >  and  if  we 

replace 


1 


with 


1 


CB  .  ^Xb,Xr. 

and  reaggregate  with  a  weighted  sum  of  squares,  then 


or  ^Bi  =  ■ 


These  are  exactly  the  ATCAL  equations! 

The  assumption  that  X b-  Xr^  is  a  time  constant  -  in  particular  that 

some  kind  of  time  to  completion  of  the  -interaction  is  a  weak  one.  For  example,  in  a 
simple  Lanchester  model: 

=  -LiiRi(<)  Ln  >  0 

ctt  \ 

=  -KuB,{t)  Ku  >  0 

at 

where  at  t  =  0,  5i(0)  >  0  and  i2i(0)  >  0,  and  at  t  =  tf,  B\(tj')  >  0  and  >  0  with 

tf  is  being  duration  of  combat,  we  would  have  the  intensity  I  =  y/KuLn  being  the 
“time  constant”  with  units  =  (time)“^. 
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Computations  yield  a  time  to  completion  of 


where 


ao  = 


bo  — 


Bijtf) 

Bi{0) 

Rijtf) 

Ri{0) 


<  1 
<  1 


are  the  fractions  remaining  at  end  of  engagement.  In  short  the  time  to  completion  is  not 
1 


In  addition,  we  recall  that  the  eigenvalue  method  depends  on  constant  kill  rates  to 
get  non- time  dependent  values.  Thus  if  our  attrition  model  is  not  of  Lan Chester  type: 


—  =  MX  M  constant 
dt 

the  values  should  either  be  time  varying  or  thought  of  as  instantaneous.  On  the  other 
hand,  ATCAL  postulates  the  attrition  as  a  function  of  time  via  its  exponential  averaging 
assumption,  given  initial  and  (derived)  final  forces.  COSAGE  does  even  better,  giving 
actual  time  dependent  loss  data.  This  allows  for  the  possibility  of  better  (time  dependent) 
values  and  (time  dependent)  shifting  of  of  fires,  and  would  move  the  ATCAL  importances 
closer  to  our  Fire  Allocation  Methods,  which  do  exhibit  shifting  allocations  of  fire  and 
changing  values  over  time.  As  it  stands,  ATCAL  gives  no  time  dependent  information; 
it  is  a  static  model  giving  the  endpoint  as  a  function  of  the  initial  point.  In  fact,  AT¬ 
CAL  postulates  an  exponential  attrition  path  between  the  initial  and  final  points,  while 
Lanchester  postulates  a  quadratic  attrition  path. 

Finally  one  can  consider  sonie  extreme  cases.  For  example,  if  we  cause  Blue  weapons 
to  kill  twice  as  many  Red  in  time  AT  then  the  ratio  of  importance  doesn’t  double.  ATCAL 
gives  y/2  as  the  relative  factor.  Another  such  case  to  consider  is  when  many  Blue  types 
oppose  one  Red  type: 

Xi,...,Xn  Blue  Xo  Red 
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n 

Set  c  =  then  the  ATCAL  importances  become: 

fc=i 


Xo=c-'/M 


3/8 


vi=l 


Xi  =  c-^/«a,o 


1/8 


diO 


.1=1 


_  <^i0  -^1/3 


Now  letting  fio  equal  the  fractional  involvement  of  X,-  in  the  kills  of  Xo,  so  that  fio 
we  get  a  relative  value 

^0 


fliQ 
C  ’ 


which  strangely  depends  on  Xq. 
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§4.  Numerical  Comparisons  and  Sensitivity  —  Results  of  Some  Simple  Numer¬ 
ical  Experiments. 

We  conducted  a  large  number  of  simple  numerical  experiments  for  each  of  the  three 
weapon  scoring  methods  discussed  above.  Our  purposes  were  two-fold.  First,  we  wanted 
to  examine  how  each  method  responded  to  various  changes  in  its  inputs.  This  led  us 
to  compute  the  amounts  of  change,  certain  numerical  derivates,  and  the  distributions 
obtained  by  treating  various  inputs  as  random  variables.  Our  objective  was  to  determine 
if  the  models  responded  in  reasonable  and  consistent  ways.  Second,  to  the  extent  possible, 
we  sought  comparisons  across  methods  that  would  allow  us  to  evaluate  the  relative  merits 
of  each  method. 

On  the  whole,  the  results  offered  no  surprises;  and  while  the  methods  produced  dif¬ 
ferent  values  for  various  weapon  systems,  the  results  were  usually  qualitatively  the  same 
for  all  three  methods.  For  that  reason,  we  shall  describe  below  only  a  small  portion  of  the 
number  of  runs  actually  done.  This  sample  is  a  good  representative  cross-section  of  the 
entire  body  of  results. 

Eigenvalue  or  Potential/Antipotential  Method. 

Following  our  notation  in  section  one,  we  take  as  initial  forces: 


and  attrition  rates: 


K  = 


L  = 


.1 

.05 

^  \ 

.075 

.1 

.025 

.025 

.025 

.075  J 

.075 

.05 

.01  \ 

.025 

.1 

.03  . 

.01 

.01 

.04  j 

Recall  that  Kij  is  the  rate  at  which  one  Blue  weapon  of  type  i  (i  =  1,2,  or  3)  kills  Red 
weapons  of  type  j  {j  =  1,2,  or  3)  and  that  Lji  is  the  rate  at  which  a  unit  of  Red  j  kills 


Blue  i. 
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The  eigenvalue  method  gives  values: 


/VB^\  /l.O  \ 

7^2  =  1.25447 

\VB3/  480896/ 

/VRi\  /. 950282  \ 
yi?2  =  1.09924 

\VRzJ  V -278556/ 


so  that  the  total  Blue  value  is  87.3208  and  the  total  Red  value  is  144.569.  This  makes  for 
a  Red  to  Blue  force  ratio  of  1.6556  to  1.  Note  that  VBi  has  been  normalized  to  be  1.  Also 
the  common  eigenvalue  A  is  .022497  and  the  intensity  /  =  y/X  =  .14999. 


Our  first  numerical  experiment  involved  computing  the  numerical  partial  derivatives 
of  the  values,  the  eigenvalue,  and  the  intensity  as  functions  of  Kij  and  Ljj.  The  results  for 
Til,  Li2,  Li3,  L21  and  T31  are  shown  below. 

deriv[K,L,NB,NR,{l,l}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 


lambda 

0.0723231 

Sqrt [lambda] :  0.241093 

VB [2]  : 

-2.4844 

VB [3]  : 

-1.13314 

VR[1]  : 

4.2359 

VR[2]  : 

-3.64993 

VR[3]  : 

-0.915579 

B:  -76.6228 
R:  259.281 
step=0 . 00001 

{0.77  Second,  Null} 

deriv[K,L,NB,NR,{l,2}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

lambda:  0.0907275 

Sqrt [lambda] :  0 . 302445 

VB[2]:  -3.11662 

VB[3]:  -1.4215 

VR[1]:  5.31383 

VR[2]:  -4.57875 

VR[3]:  -1.14857 

B:  -96.1214 

R:  325.261 

step=0 .00001 

{0.66  Second,  Null) 
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deriv[K,L,NB,NR,{l,3}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

lambda:  0.0347799 

Sqrt [lambda] :  0.115941 

VB[2]:  -1.19474 

VBC3]:  -0.544925 

VR[1]:  2.03703 

VR[2]:  -1.75524 

VRC3]:  -0.440299 

B:  -36.8476 

R:  124.687_ 

step=l.  10”° 

{0.77  Second,  Null} 

deriv[K,L,NB,NR,{2,l}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

lambda:  0.0677829 

Sqrt [lambda] :  0.225958 

VB[2]:  1.77313 

VB[3]:  0.0206782 

VR[1]:  -0.839129 

VR[2]:  6.19742 

VR[3]:  -0.295911 

B:  42.8652 

R:  158.066 

step=0 . 00001 

{0.71  Second,  Null} 

deriv[K,L,NB,NR,{3,l}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

lambda:  0.0242435 

Sqrt [lambda] :  0.0808172 

VB[2]:  1.47834 

VB[3]:  3.78407 

VR[1]:  0.233072 

VR[2] :  1.1502 

VR[3]:  7.62474 

B:  92.2411 

R:  221.81 

step=l.  10~® 

{0.65  Second,  Null} 


For  example  increasing  L12  amounts  to  making  Red  1  more  effective  against  Blue 
2.  As  a  result  the  value  of  Red  1,  VR\^  increases,  as  does  the  eigenvalue  A,  the  intensity 
I  =  \/A,  and  the  total  Red  value.  All  other  values,  including  the  total  Blue  value,  decrease. 
At  first  this  might  seem  odd,  but  the  reader  is  urged  to  keep  in  mind  that  these  values 
are  relative  values  (relative  to  JBj),  not  absolute  values,  and  this  complicates  the  situation. 
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Recall  also  that  L12  =  -05  so  that  a  10%  increase  in  effectiveness  of  Red  1  against  Blue  2 
amounts  to  an  increase  of  .005  in  L12.  This  decreases  the  value  of  a  unit  of  Blue  2  relative 
to  a  unit  of  Blue  1  by  3.11662  times  .005  or  .0155831.  This  is  a  relatively  small  change. 

Next,  we  did  a  direct  check  on  the  numerical  derivative  by  altering  (increasing)  the 
Kij  and  Lji  entries  by  1%,  5%  and  10%  and  recomputing  all  quantities.  Sample  results 


appear  below; 

compentry[K,L,NB,NR,{l,l5}] 


Actual 

+1% 

+5% 

+10% 

lambda 

0.022497 

0.022551 

0.02277 

0.023045 

VB[2] 

1.2545 

1.2526 

1.2453 

1.2363 

VB[3] 

0.4809 

0.48005 

0.4767 

0.47262 

VR[1] 

0.95028 

0.95346 

0.9661 

0.9818 

VR[2] 

1.0992 

1.0965 

1.0857 

1.0725 

VR[3] 

0.27856 

0.27787 

0.27516 

0.27185 

B 

87.321 

87.264 

87.037 

86.761 

R 

144.57 

144.76 

145.54 

146.52 

compentry[K,L,NB,NR,{l,2}] 

Actual 

+1% 

+5% 

+10% 

lambda 

0.022497 

0.022542 

0.022723 

0.022949 

VB[2] 

1.2545 

1.2529 

1.2468 

1.2394 

VB[3] 

0.4809 

0.48019 

0.4774 

0.47403 

VR[1] 

0.95028 

0.95293 

0.96344 

0.97636 

VR[2] 

1.0992 

1.097 

1.088 

1.0771 

VR[3] 

0.27856 

0.27798 

0.27573 

0.27299 

B 

87.321 

87.273 

87.085 

86.856 

R 

144.57 

144.73  ^ 

145.38 

146.18 
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compentry[K,L,NB,NR,{2,2}] 


Actual 

+1% 

+5% 

+10% 

lambda 

0.022497 

0.022582 

0.022924 

0.023353 

VB[2] 

1.2545 

1.2567 

1.2655 

1.2762 

VB[3] 

0.4809 

0.48092 

0.48103 

0.48116 

VR[1] 

0.95028 

0.94923 

0.94504 

0.93983 

VR[2] 

1.0992 

1.107 

1.138 

1.1767 

VR[3] 

0.27856 

0.27819 

0.27671 

0.27489 

B 

87.321 

87.374 

87.587 

87.847 

R 

144.57 

144.77 

145.56 

146.55 

comprow[K,L,NB,NR,l] 

Actual 

+1% 

+5% 

+10% 

lambda 

0.022497 

0.0226 

0.023014 

0.023534 

VB[2] 

1.2545 

1.251 

1.2373 

1.2212 

VB[3] 

0.4809 

0.47929 

0.47307 

0.46572 

VR[1] 

0.95028 

0.9563 

0.98005 

1.009 

VR[2] 

1.0992 

1.0941 

1.074 

1.0501 

VR[3] 

0.27856 

0.27726 

0.27222 

0.26623 

B 

87.321 

87.212 

86.792 

86.295 

R 

144.57 

144.94 

146.41 

148.23 

compcol[K,L,NB,NR. 

.1] 

Actual 

+1% 

+5% 

+10% 

lambda 

0.022497 

0.022571 

0.022866 

0.023238 

VB[2] 

1.2545 

1.2532 

1.2482 

1.2422 

VB[3] 

0.4809 

0.48043 

0.47861 

0.4764 

VR[1] 

0.95028 

0.95327' 

0.96516 

0.97989 

VR[2] 

1.0992 

1.0982 

1.094 

1.089 

VR[3] 

0.27856 

0.27856 

0.27859 

0.27865 

B 

87.321 

87.283 

87.137 

86.959 

R 

144.57 

144.82 

145.85 

147.12 
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Note  that  comprow[K,L,NB,NR,l]  changes  row  1  of  L  increasing  in,  in,  ii3  all  by 
the  same  percentage.  This  corresponds  to  making  Red  1  more  effective  against  all  Blue 
types.  Similarly  compcol[K,L,NB,NR,l]  changes  column  1  of  i  increasing  in,  i2i,  isii 
in  effect  making  all  Red  weapons  more  effective  against  Blue  1. 

Finally,  we  treated  various  entries  and  Lji  as  random  variables.  Besides  varying 
them  one  at  a  time,  we  varied  whole  columns  and  rows  of  K  and  L  by  varying  the  entries 
both  independently  and  together  (correlated);  We  also  treated  the  entire  matrix  A  and 
i  as  random.  We  used  both  normal  distributions  and  beta  distributions  to  model  the 
stochastic  behavior. 

Typical  results  are  shown  below.  In  this  first  example  Lii  is  treated  as  a  random 
variable  subject  to  a  normal  distribution  with  mean  .075  (the  original  value  of  An)  and 
standard  deviation  equal  to  5%  of  the  mean  value,  i.e.  .00375.  The  program  did  1500 
draws  from  this  distribution  and  calculated  the  distributions  for  the  Red  and  Blue  values. 

normentry[K,L,NB,NR,1500,.05,{l,l}]/ /Timing 

Lambda:  0.022497 

Mean:  0.0225001 

Median:  0.0225026 

Variauice:  7.07888  10 

Standard  Deviation:  0.000266062 

REinge:  0.00163993 

Minimum:  0.0216807 

Maximum:  0.0233206 

Sqrt [Lambda] :  0 . 14999 
Mean:  0.149998 
Median:  0.150009  _ 

Variance:  7.86735  10” 

Stamdard  Deviation:  0.000886981 
Range:  0.00546728 
Minimum :  0 . 147244 
Maximum:  0.152711 

VB[2]:  1.25447 
Mean:  1.25453 
Median:  1.25428 
Variance:  0.0000837406 
Standard  Deviation:  0.00915099 
Range:  0.0564821 
Minimum:  1.22768 
Maximum :  1 . 28417 
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VB[3]:  0.480896 

Mean:  0.480924 

Median:  0.480809 

Variance:  0.0000174203 

Standard  Deviation:  0.00417377 

Range:  0.0257613 

Minimum :  0 . 468676 

Maximum :  0 . 494437 

VR[1]:  1. 

Mean:  0.950323 

Median:  0.950607 

Variance:  0.000243054 

Standard  Deviation:  0.0155902 

Range:  0.0961508 

Minimum:  0.901104 

Maximum:  0.997255 

VR[2]:  1.25447 
Mean:  1.09931 
Median:  1.09896 
Variaoice:  0.00018069 
Standard  Deviation:  0.0134421 
Range:  0.0829561 
Minimum:  1.05971 
Maximum :  1 . 14266 

VRC3]:  0.480896 

Mean:  0.278573 

Median:  0.278486 

Variaoice:  0.0000113697 

Standard  Deviation:  0.0033719 

Range:  0.020809 

Minimum :  0 . 268638 

Maximum :  0 . 289447 

B:  87.3208 

Mean:  87.3227 

Mediein:  87.3149 

Variauice:  0.0796539 

Standard  Deviation:  0.28223 

Range:  1.74199 

Minimum :  86 . 4946 

Maximum :  88 . 2366 

R:  144.569 

Meeoi:  144.576 

Mediaui:  144.589 

Variance:  0.910026 

Standard  Deviation:  0.953953 

Range:  5.88066 

Minimum:  141.606 

Majcimum :  147 . 487 


VR[2]  Distribution 


86.336.75  87  87.2387.337.75  88 
R  Distribution 


142  143  144  145  146  147 


In  order  to  make  comparisons  with  ATCAL,  it  was  necessary  to  create  examples 
where  the  killer- victim  scoreboard  yielded  essentially  the  same  kill  rates  as  those  used  in 
the  eigenvalue  examples.  To  some  extent  this  is  a  guess,  in  that  the  attrition  model  is 
unspecified.  We  simply  assumed  straight  line  attrition. 

For  example,  take  initial  forces  of 

/BA  /2000\  /  Ri\  /4000\ 

B2  I  =  960  and  i22  =  1600 

ynj  \60oJ  \R3/  \80oJ 


with  killer  victim  scoreboard 


killer 


victim 


Bx  B2  Bz  Ri  R2  Rz 

Bi  0  0  0  200  100  0 

B2  0  0  0  72  96  24 

Bz  0  0  0  15  15  45 

Ri  300  200  40  0  0  0 

R2  40  160  48  0  0  0 

R3  8  8  32  0  0  0 


Assuming  a  duration  of  1  time  unit  and  straight  line  attrition,  we  get 


K  = 


/.I  .05  0  \  /.075  .05  .01  \ 

.075  .1  .025  L  =  [  .025  .01  .03 

\.025  .025  .075/  V-01  / 


as  we  had  in  the  previous  eigenvalue  method  examples.  Using  ATCAL,  we  obtain  relative 
importances  of 

VBi  =1  VRi=  .94595 

VB2  =  1.5513  VR2  =  1.5396 

VB3  =  1.7971  VRz  =  1.0477 

This  result  appears  to  differ  from  the  eigenvalue  result  in  that  both  Red  and  Blue  type  3 
appear  more  valuable. 

The  numerical  derivatives  with  respect  to  varying  various  Lji  (or  equivalently  vary¬ 
ing  entries  in  the  killer-victim  scoreboard)  are  shown  below.  We  give  both  the  absolute 
derivatives  and  the  relative  derivatives  with  VBi  normalized  to  1. 


67 


F 


BK={{200,100,0}, {72,96,24},  115,15,45}}; 
RK={{300,200,40},  {40,160,48},  {8,8,32}}; 
NB={2000,960,600}; 
NR={4000,1600,800}; 


deriv[BK,RK,NB,NR,{l,l}//Timing 

Approximate  Numerical  Derivatives  (5  point) 
VB[1]:  0.0000373978 
VB[2]:  6.85268  10“® 

VB[3]:  1.31578  lO""^ 

VRCl]:  0.000179095 
VR[2]:  1.07447  10"® 

VR[3]:  -9.29081  10"® 

{9.5  Second,  Null} 

deriv[BK,RK,NB,NR,{l,2}]//Timing 


Approximate  Numerical  Derivatives  (5  point) 


VB[1] 

VB[2] 

VB[3] 

VR[1] 

VR[2] 

VR[3] 


0.0000393462 
-0.0000469537 
-6.93696  10"”^ 
0.000270991 
-0.000193113 
-5.39146  10"'^ 


{9.23  Second,  Null} 


deriv[BK,RK,NB,NR,{l,3,}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

VB[1]:  3.38649  10"® 

VBC2]:  -0.0000430005 
VBC3]:  -0.000275548 
VR[1]:  0.0000618685 
VR[2]:  -0.000107301 
VR[3]:  -0.00048685 

{9.45  Second,  Null} 

deriv[BK,RK,NB,NR{2,l}]//Timing 


Approximate  Numerical  Derivatives 
VB[1]:  -5.0077  10 


-6 


VBC2] 
VB  [3] 
VR[1] 
VR[2] 
VRC3] 


1.99616  10"® 
-2.68284  10"® 
-0.0000284272 
0.0000103152 
-5.92568  10"® 


{11.09  Second,  Null} 


(5  point) 
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deriv[BK,RK,NB,NR,{3,l}]//Timing 
Approximate  Numerical  Derivatives  (5  point) 
VB[1]:  -3.02879  10"® 

VB[2];  -4.93462  10"^ 

VB[3]:  2.19239  10“® 

VR[1] :  -0.0000144314 
VR[2]:  -2.592  lO""^ 

VR[3];  3.90199  10"® 

{13.34  Second,  Null} 
deriv[BK,RK,NB,NR,{2,2}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 
VB[1]:  0.000040094 
VB[2]:  0.000136194 
VB[3]:  1.74046  10"® 

VR[1]:  -4.4331  10"® 

VR[2]:  0.000463299 
VR[3]:  2.78241  lO""^ 

(12.42  Second,  Null} 

BK={{200,100,0},  {72,96,24},  (15,15,45}}; 
RK={{300,200,40},  {40,160,48},  {8,8,32}}; 
NB={2000,960,600}; 

NR={4000,1600,800}; 

deriv[BK,RK,NB,NR,{l,l,}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 


VB[2]  : 

-1.37032 

VB [3]  : 

-1.79654 

VR[1]  : 

3.84922 

VR[2]  : 

-1.51334 

VR[3]  : 

-1.05191 

{12.14  Second,  Null} 
deriv[BK,RK,NB,NR,{l,2}]//Timing 
Approximate  Numerical  Derivatives  (5  point) 


VB  [2]  ; 

-2.89237 

VB [3]  : 

-1.91243 

VRCl]  : 

6.2611 

VR[2]  : 

-6.79461 

VR[3]  : 

-1.11854 

{13.78  Second,  Null} 
deriv[BK,RK,NB,NR,{l,3}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

VB[2]:  -1.29239 

VB[3]:  -7.54303 

VR[1]:  1.57123 

VR[2]:  -3.0135 

VRC3]:  -13.1236 

{20.55  Second,  Null} 
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deriv[BK,RK,NB,NR,{2,l}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

VB[2]:  0.26153 

VB[3];  0.240317 

VR[1] :  -0.634494 

VR[2]:  0.482767 

VRC3]:  0.138934 

{9.45  Second,  Null} 

deriv[BK,RK,NB,NR,{3,l}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

VBC2]:  0.112628 

VBC3]:  0.204504 

VR[1];  -0.309782 

VRC2]:  0.117952 

VR[3];  0.189498 

{9.44  Second,  Null) 

deriv[BK,RK,NB,NR{2,2}]//Timing 

Approximate  Numerical  Derivatives  (5  point) 

VBC23:  1.98182 

VB[3]:  -1.88323 

VRCl]:  -1.13453 

VR[2];  10.7553 

VRC3]:  -1.11763 

{9.4  Second,  Null) 


For  the  Fire  Allocation  method,  we  follow  the  simple  model  as  outlined  in  Sec 
tion  2.a.ii.  NB  and  NR  represent  the  initial  forces  and  A  and  C  are  the  pure  kill  rates 
following  the  notation  of  Section  2.a.ii,  except  that  the  indices  on  aij  and  Cij  have  been 
switched.  Thus  021  in  the  program  is  012  in  the  text  of  Section  2.a.ii.  Note  also  that  the 
value  hi  is  normalized  to  1. 
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Allocation  Method 
Base  Program 
June  3, 1994 

This  program  determines  the  weapon  values  given  the  matrices  A  and  C  and 
the  strength  of  the  forces  NB  and  NR. 

Off [Syntax; : newl , General : : spell , General : : spelll , Part: :partd] 

allocate [ { { all_, al2_} , { a21_,  a22_}  } , { { cll_,  cl2 _ } , {c21 _ ,  c22_}  } , 

{b01_,b02_},{r01_,r02_}]  :  = 

Block [ { listl , list2 , lists , list4 ,phi ,psi ,matx,philist ,psilist , tl , 
t2 , L  f  j , common, lam, X, list ^ xvec) , 

listl={-all  b01''2  cll''2  r01+b01^2  cll  phi  r01-al2  b01''2  c21''2  r02 
+b01^2  c21  phi  r02,  bOl  b02  cll  phi  rOl  +  bOl.  b02  c21  phi 

r02,  -2  all  bOl  b02  cll  cl2  rOl  +  bOl  b02  cl2  phi  rOl  -  2  al2 

bOl  b02  c21  c22  r02  +  bOl  b02  c22  phi  r02,  b02^2  cl2  phi  rOl 
+  b02^2  c22  phi  r02,  -all  b02"'2  cl2^2  rOl  -  al2  b02^2  c22^2 
r02 , 0 , 0 , 0 } ; 
matx={ } ; 

For  [i=0,i<=3,i++,i^pendTo[inatx,RotateRight[listl,i]  ]  ]  ; 

Iist2={-a21  b01-'2  cll^2  rOl  -  a22  b01^2  c21^2  r02, 

b01^2  cll  phi  rOl  +  b01^2  c21  phi  r02,  -2  a21  bOl  b02  cll  cl2 

rOl  +  bOl  b02  cll  phi  rOl  -  2  a22  bOl  b02  c21  c22  r02  +  bOl 

b02  c21  phi  r02,  bOl  b02  cl2  phi  rOl  +  bOl  b02  c22  phi  r02, 
-a21  b02-^2  cl2^2  rOl  +  b02^2  cl2  phi  rOl  -  a22  b02''2  c22^2 
r02  +  b02^2  c22  phi  r02, 0,0,0} ; 

For [i=0 , i<=3 , i++ ,2^pendTo [matx, RotateRight [list2 , i ] ] ] ; 
philist=Select  [phi  /.  NSolve [Det  [matx]  =0, phi]  ,  Im[#]—0  6]  ; 
list3={ -ally's  bOl  cll  r01^2  -  a21^2  b02  cl2  r01^2  +  all  bOl 
psi  r01''2  +  a21  b02  psi  r01''2,  all  bOl  psi  rOl  r02  + 
a21  b02  psi  rOl  r02,  -2  all  al2  bOl  cll  rOl  r02  -  2  a21  a22 

b02  cl2  rOl  r02  +  al2  bOl  psi  rOl  r02  +  a22  b02  psi  rOl  r02, 

al2  bOl  psi  r02''2  +  a22  b02  psi  r02^2,  -al2^2  bOl  cll  r02^2  - 
a22^2  b02  cl2  r02^2,0,0,0} ; 
inatx={  } ; 

For [i=0,i<=3,i++,  i^pendTo [matx, RotateRight [lists, i] ] ] ; 

list4={-all''2  bOl  c21  rOl^^S  -  a21'^2  b02  c22  r01^2,  all  bOl  psi 

r01'‘2  +  a21  b02  psi  r01^2,  -2  all  al2  bOl  c21  rOl  r02  -  2  a21 

a22  b02  c22  rOl  r02  +  all  bOl  psi  rOl  r02  +  a21  b02  psi  rOl 

r02,  al2  bOl  psi  rOl  r02  +  a22  b02  psi  rOl  r02,  -al2^2  bOl 
c21  r02''2  -  a22''2  b02  c22  r02''2  +  al2  bOl  psi  r02''2  +  a22  b02 
psi  r02''2,0,0,0} ; 

For [i=0 , i<=3 , i++ ,  i^pendTo [matx, RotateRight [list4 , i ] ] ] ; 
psilist=Select  [psi  /.  NSolve  [Det  [matx]  **0 ,  psi] ,  Im[#]—0  4], 

coinmon={ } ; 
tl=Length[philist] ; 
t2=Length[psilist] ; 

For[i=l,i<=tl,i++, 
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For [j=l, j<=t2, j++, 

If[  Abs[philist[ [i] ]-psilist[ [j] ]  ]<10^-10, 

i^pendTo [ common , philis t [ [ i ] ]  ] 

]; 

]  ; 

]  ; 

If  [ Length  [  common ]  =0 ,  Break  [  ]  ]  ; 

lam=Max[ common] ; 

Print [ "Common  Lambda  Value:  ",1am]; 

xvec={x''4 ,  x''3 ,  x''2  ,  x,  1 }  ; 
list=Take[listl,5]  /.  phi->lam; 

xvalsl=Select  [x  /.  Solve  [list.  xvec=0,  x]  ,  lm[#]=0  &]; 

list=Take[list2,5]  /.  phi->lam; 

xvals2=Select  [x  /.  Solve  [list .  xvec=0  ,x]  ,  lm[#]=0  &]  ; 
common* { } ; 
tl=Length [xvalsl] ; 
t2=Length[xvals2] ; 

For [i=l , i<=tl , i++ , 

For [ j=l , j<=t2 , j++ , 

If[  Abs  [xvalsl  [  [i]  ]-xvals2[  [j]  ]  ]<10''-10, 

AppendTo [ common , xval s 1 [ [ i ] ]  ] 

]  ; 

]  ; 

]  ; 

If  [  Length  [  common  ]  =0 ,  Break  [  ]  ]; 

ratio=Max[ common] ; 

hl=l  ; 

h2=l /ratio; 

gl=(b01  cll  hl''2  +  b02  cl2  h2''2)  /  (Sqrt  [lam]  *  (bOl  hi  +  b02  h2) ) 
g2=(b01  c21  hl''2  +  b02  c22  h2''2)  /  (Sqrt  [lam]  *  (bOl  hi  +  b02  h2) ) 

Print["H[l]=",hl] ; 

Print ["H [2]=" ,h2] ; 

Print["G[l]*" ,gl] ; 

Print["G[2]=" ,g2] 

] 

allocate [{{ .1, .05},{ .075, .!}},{{ .075, .05}, { . 025 ,.!}},{ 50 , 25 } , 
{100,40}] 

Common  Lambda  Value:  0.00561786 
H[l]=l 

H[2]=0. 849138 
G[1]=0. 871234 
G[2]-0. 571781 


Allocation  Method 
Numerical  Derivatives 
June  3, 1994 

This  program  accapts  the  the  input  data  of  the  A  and  C  matrices,  the  strength  of 
the  forces  NB  and  NR,  and  the  indices  of  an  element  in  the  C  matrix,  returning  the 
numerical  derivatives  of  the  weapon  values  with  respect  to  this  element. 

off [Syntax: : newl , General : : spell , General : : spelll , Part : :partd] 
deriv[A_,C_,NB_,NR_,ind_]  :  = 

Block [ [ NA, NC , step , Rlist , Blist , base , process , bOl , b02 , rOl , r02 , x, 
xvec, vals} , 

base  [inatA_,iiiatC_]  :  = 

Block [ [ all , al2 , a21 , a22 , cll , cl2 , c21 , c22 , listl , list2 , lists , list4 , 
matx,  i ,  j  »phi  ,philist ,psi  ,psilist,  common ,  tl,,  t2  ,  lam,  list, 
xvalsl , xvals2 , ratio , gl , g2 , h2 } , 

{  {all,al2}  ,  {a21,a22}  }=IaatA; 

{{cll,cl2},{c21,c22}}=matC; 

Inatx={  } ; 

listl={-all  b01'*'2  cll''2  r01+b01''2  cll  phi  r01-al2  b01^2  c21‘'2  r02 
+b01''2  c21  phi  r02,  bOl  b02  cll  phi  rOl  +  bOl  b02  c21  phi 

r02 ,  -2  all  bOl  b02  cll  cl2  rOl  +  bOl  b02  cl2  phi  rOl  -  2  al2 

bOl  b02  c21  c22  r02  +  bOl  b02  c22  phi  r02,  b02^2  cl2  phi  rOl 
+  b02^2  c22  phi  r02,  -all  b02''2  cl2''2  rOl  -  al2  b02''2  c22"'2 
r02,0,0,0>; 

For  [i=0,i<=3,i++, AppendTo[matx, Rota teRight[listl,i] ] ] ; 

Iist2={-a21  b01''2  cll^2  rOl  -  a22  b01''2  c21''2  r02 , 

b01''2  cll  phi  rOl  +  b01''2  c21  phi  r02,  -2  a21  bOl  b02  cll  cl2 

rOl  +  bOl  b02  cll  phi  rOl  -  2  a22  bOl  b02  c21  c22  r02  +  bOl 

b02  c21  phi  r02,  bOl  b02  cl2  phi  rOl  +  bOl  b02  c22  phi  r02 , 
-a21  b02''2  cl2''2  rOl  +  b02''2  cl2  phi  rOl  -  a22  b02^2  c22''2 
r02  +  b02^2  c22  phi  r02, 0,0,0 }; 

For [ i=0 , i<=3 , i++ ,  ;^pendTo [matx , RotateRight [ lis t2 , i ] ] ] ; 
philist=Select  [phi  /.  NSolve  [Det  [matx]=0,phi]  ,  lm[#]=0  fi]  ; 


matx={ } ; 

list3={-all^2  bOl  cll  r01^2  -  a21''2  b02  cl2  r01^2  +  all  bOl 
psi  r01^2  +  a21  b02  psi  r01^2,  all  bOl  psi  rOl  r02  + 
a21  b02  psi  rOl  r02,  -2  all  al2  bOl  cll  rOl  r02  —  2  a21  a22 

b02  cl2  rOl  r02  +  al2  bOl  psi  rOl  r02  +  a22  b02  psi  rOl  r02, 

al2  bOl  psi  r02''2  +  a22  b02  psi  r02^2,  -al2^2  bOl  cll  r02^2  - 
a22''2  b02  cl2  r02''2', 0 , 0 , 0 }  ; 

For  [i=0 , i<=3 , i++ , AppendTo [matx, RotateRight [lists , i ] ] ] ; 
list4={-all^2  bOl  c21  rOl^^S  -  a21''2  b02  c22  rOl-'S,  all  bOl  psi 

r01'‘2  +  a21  b02  psi  r01^2,  -2  all  al2  bOl  c21  rOl  r02  -  2  a21 

a22  b02  c22  rOl  r02  +  all  bOl  psi  rOl  r02  +  a21  b02  psi  rOl 

r02,  al2  bOl  psi  rOl  r02  +  a22  b02  psi  rOl  r02 ,  -al2''2  bOl 
c21  r02^2  -  a22''2  b02  c22  r02''2  +  al2  bOl  psi  r02'*‘2  +  a22  b02 
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psi  r02''2 , 0 , 0 , 0  }  ; 

For [i=0,i<=3,i++,i^pendTo[matx,RotateRight[list4,i] ] ] ; 
psilist=Select  [psi  /.  NSolve[Det  [niatx]==0,psi]  ,  lm[#]==0 

conimon=  {  }  ; 
tl=Length[philist] ; 
t2=Length [psilist] ; 

For [i=l , i<=tl , i++ , 

For[j=l, j<=t2, j++, 

If[  Abs[philist[ [i] ]-psilist [ [ j] ]  ]<10^-10, 

i^pendTo [common, philist [ [i] ]  ] 

]  ; 

]  ; 

]  ; 

If  [Length  [  common] —0 , 

Print ["***  NO  COMMON  LAMBDA  VALUE?  ***"]; 

Break [ ] 

]  ; 

If [Length [ common] >1 , 

Print  [ "  *  *  *  MULTIPLE  COMMON  LAMHDA  VALUES !  *  *  *  "  ] 

]  ; 

lam=Max[ common] ; 

list=Take[listl,5]  /.  phi->lam; 

xvalsl=Select[x  /.  Solve[list.xvec~0,x]  ,lm[#]==0  S]; 
list=Take[list2,5]  /.  phi->lam; 

xvals2*Select [x  /.  Solve[list.xvec=0,x]  ,lm[#]=0  &]  ; 
common= { }  ; 
tl=Length[xvalsl] ; 
t2=Length[xvals2] ; 

For [i=l , i<=tl , i++ , 

For [ j=l , j<=t2 , j++ , 

If[  Abs  [xvalsl  [  [i]  ] -xvals2  [  [  j  ]  ]  ]<10''-10, 

j^pendTo [ common , xvalsl [ [ i ] ]  ] 

]  ; 

]  ; 

]; 

If  [Length  [ common]  —0 , 

Print["***  NO  ccadMON  RATIO?  ***'']; 

Break [ ] 

]; 

If [Length [ common] >1 , 

Print ["***  MULTIPLE  COMMON  RATIOS?  ***"] 

]; 

ratio=McCC  [ common]  ; 
h2=l /ratio; 

gl=(b01  cll  +  b02  cl2  h2''2)  /  (Sqrt  [lam]  *  (bOl  +  b02  h2) ) 
g2=(b01  c21  +  b02  c22  h2^2) / (Sqrt [lam] * (bOl  +  b02  h2) ) 
{{gl,g2},h2} 


]  ; 

process [vals_]  :  = 

Block [ { } , 

AppendTo [Blist, vals [ [2] ]  ]; 

i^pendTo  [Rlist  [  [1]  ]  ,vals  [[1,1]]  ]  ,’ 

AppendTo [Rlist [ [2 ]], vals [ [1 , 2 ] ]  ]; 

]  ; 

calc[x_] :=(x[[l]]-8x[[2]]+8x[[3]]-x[[4]])/(12  step) ; 

NA=A; 

NC=C ; 
step=l ; 

While[NC[[ind[[l]] ,ind[[2]]  ] ] -2*step<=0 , step*=. 1] ; 

step*=.001; 

{b01,b02}=NB; 

{r01,r02}=NR; 
xvec={x''4  ,x‘'3,x''2  ,x,  1}  ; 

Rlist={ {},{}}; 

Blist={ } ; 

NC[ [ind[ [1] ] ,ind[ [2]]  ]]-=2*step; 

process [base [NA,NC] ] ; 

NC[ [ind[ [1] ] ,ind[ [2] ]  ]]+=step; 
process [base [NA,NC] ] ; 

NC[ [ind[ [1] ] ,ind[ [2] ]  ] ]+=2*step; 

process [base [NA,NC] ] ; 

NC[ [ind[ [1] ] ,ind[ [2] ]  ]]+=step; 

process [base [NA,NC] ] ; 

Print ["^^proximate  Numerical  Derivatives"]; 

Print["G[l]:  ",  calc [Rlist [ [1 ]] ]  ]; 

Print ["G [2]:  ",  calc [Rlist [ [2 ]] ]  ]; 

Print["H[2]:  " , calc[Blist]  ] ; 

] 

A={{.1,.05},{.075,.1}}; 

Cx={{ .075, .05} , { .025, .1}}; 

NB={ 50,25 } ; 

NR={100,40}; 

Derivatives  with  respect  to  C[1,1] 

deriv[A,Cx,NB,NR,  {1,1}] 

Approximate  Numerical  Derivatives 
G[l] :  4.29571 
G[2] :  -4.59138 
H[2] :  -2.36597 

Derivatives  with  respect  to  C[1,2] 

deriv[A,Cx,NB,NR, {1,2}] 

Approximate  Numerical  Derivatives 
G[l] :  1.54868 
G[2] :  -1.65528 

H[2]  ;  -0.852976 _ _ 

Numdcriv  3 
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Derivatives  with  respect  to  C[2,1] 

deriv [A, Cx,NB,NR,  {2,1}] 

Approximate  Numerical  Derivatives 
G[l]:  0.589737 
G[2] :  11.6787 
H[2] ;  3.60508 

Derivatives  with  respect  to  C[2,2] 

deriv  [A,  Cx,NB,NR,  {2,2}] 

Approximate  Numerical  Derivatives 
G[l] :  0.212611 
G[2]  :  4.2104 
H[2]:  1.2997 
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Allocation  Method 
Entry  Sensitivity  Analysis 
May  27, 1994 

This  program  receives  the  matrices  A  and  C,  the  strength  of  the  forces  NB  and  NR, 
the  number  of  iterations,  the  standard  deviation  of  the  distribution,  and  the  indices 
of  the  entry  in  the  C  matrix  to  be  altered.  First,  the  program  calculates  the  weapon 
values  for  the  initial  A  and  C.  The  program  then  repeatedly  perturbs  the  indicated 
entry  by  a  random  draw  from  the  given  normal  distribution.  The  statistics  of  the 
calculated  weapon  values  are  returned  and  the  distributions  are  plotted. 

off [General : : spelll ,  Gteneral : : spell , Syntax: : newl , Part: :partd] 
«Statistics '  DescriptiveStatistics ' 

«Statistics '  ContinuousDistributions ' 

normentry  [A_,  C_,NB_,NR_,  iter_,  stdev_,  ind_]  :  = 

Block [ { NA, NC , Blis t , Rlist , B , R, bO 1 , b02 , rO 1 , r02 , base , process , x , xvec , 
ndist  f results , i } / 

base  [matA__,matC_]  :  = 

Block [ { all , al2 , a21 , a22 , cll , cl2 , c21 , c22 , listl , list2 , list3 , list4 , 
matx, i ; j  f pbi , philist , psi , psilist , common , tl , t2 , lam, list , 
xvalsl , xvals2 , ratio , gl , g2 , h2 } , 

{ {all,al2} , {a21,a22} }=matA; 

{{cll,cl2},{c21,c22})=matC; 

matx={ } ; 

listl={-all  b01^2  cll^2  r01+b01''2  cll  phi  r01-al2  b01''2  c21''2  r02 
+b01''2  c21  phi  r02,  bOl  b02  cll  phi  rOl  +  bOl  b02  c21  phi 

r02,  -2  all  bOl  b02  cll  cl2  rOl  +  bOl  b02  cl2  phi  rOl  -  2  al2 

bOl  b02  c21  c22  r02  +  bOl  b02  c22  phi  r02,  b02''2  cl2  phi  rOl 
+  b02^2  c22  phi  r02,  -all  b02''2  cl2''2  rOl  -  al2  b02^2  022^*2 
r02 ,0,0,0}; 

For [i=0,i<=3,i++, AppendTo [matx, Rota teRight[ listl, i] ] ] ; 

Iist2={-a21  b01''2  cll''2  rOl  -  a22  b01''2  c21''2  r02, 

b01''2  cll  phi  rOl  +  b01''2  c21  phi  r02,  -2  a21  bOl  b02  cll  cl2 

rOl  +  bOl  b02  cll  phi  rOl  -  2  a22  bOl  b02  c21  c22  r02  +  bOl 

b02  c21  phi  r02,  bOl  b02  cl2  phi  rOl  +  bOl  b02  c22  phi  r02, 
-a21  b02^2  cl2^2  rOl  +  b02''2  cl2  phi  rOl  -  a22  b02''2  c22''2 
r02  +  b02^2  c22  phi  r02, 0,0,0} ; 

For [i=0 , i<=3 , i++ , AppendTo [matx ,  Rota teRight [list2 , i ] ] 1 ; 
philist=Select [phi  /.  NSolve[Det[matx]— 0,phi]  ,lm[#]=*0  fi] ; 

matx={ } ; 

list3={-all^2  bOl  cll  r01''2  -  a21''2  b02  cl2  r01^2  +  all  bOl 
psi  r01''2  +  a21  b02  psi  r01^2,  all  bOl  psi  rOl  r02  + 
a21  b02  psi  rOl  r02,  -2  all  al2  bOl  cll  rOl  r02  -  2  a21  a22 
b02  cl2  rOl  r02  +  al2  bOl  psi  rOl  r02  +  a22  b02  psi  rOl  r02, 
al2  bOl  psi  r02''2  +  a22  b02  psi  r02''2,  -al2''2  bOl  cll  r02''2  - 
a22''2  b02  cl2  r02''2,0,0,0}; 

For [i=0,i<=3,i++,i^pendTo [ma tx, Rota teRight [list3,i] ] ] ; 
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Iist4={-all''2  bOl  c21  r01^2  -  a21^2  b02  c22  r01^2,  all  bOl  psi 

r01^2  +  a21  b02  psi  r01‘'2,  -2  all  al2  bOl  c21  rOl  r02  -  2  a21 
a22  b02  c22  rOl  r02  +  all  bOl  psi  rOl  r02  +  a21  b02  psi  rOl 
r02 ,  al2  bOl  psi  rOl  r02  +  a22  b02  psi  rOl  r02,  -al2''2  bOl 
c21  r02''2  -  a22''2  b02  c22  r02'*‘2  +  al2  bOl  psi  r02^2  +  a22  b02 
psi  r02''2 , 0 , 0, 0}  ; 

For[i=0,i<=3,i++,AppendTo[matx,RotateRight[list4,i] ] ] ; 

psilist=Select [psi  /.  NSolve [Det [matx]— 0 ,psi]  , lin[#]=0  &]  ; 

coinmon=  [ } ; 
tl=Length [philist] ; 
t2=Length[psilist] ; 

For [i=l , i<=tl , i++ , 

For [ j=l , j<=t2 , j++ , 

If[  Abs[philist[  [i]  ]-psilist[  [j]  ]  ]<10''-10, 

i^pendTo  [common, philist: [  [i]  ]  ] 

]  ; 

]; 

]  ; 

If  [Length  [ common]  =0 , 

Print["***  NO  COMMON  LAMBDA  VALUE?  **♦■']; 

Break [ ] 

]; 

If [Length [ common] >1 , 

Print ["***  MULTIPLE  COMMON  LAMBDA  VALUES?  ***"] 

]  ; 

lam=MeLX [common]  ; 

list=Take [listl , 5]  /.  phi->lam; 

xvalsl=Select  [x  /.  Solve  [list .  xvec=0  ,x]  ,  lm[#]==0  &]  ; 
list=TaJce[list2,5]  /.  phi->lam; 

xvals2=Select  [x  /.  Solve  [list .  xvec=0  ,x]  ,  lm[#]=0  &]  ; 
common= { } ; 
tl=Length[xvalsl] ; 
t2=Length[xvals2] ; 

For [ i=l , i<=tl , i++ , 

For [ j=l , j<=t2 , j++ , 

If[  Abs[xvalsl[  [i]  ]-xvals2  [  [  j]  ]  ]<10''-10, 

i^pendTo  [common, xv^Q. si  [  [i]  ]  ] 

1; 

]; 

]  ; 

If  [Length  [ common]  —0 , 

Print[’'***  NO  COMMON  RATIO?  ***"]; 

Break [ ] 

]  ; 

If [Length [ common] >1 , 

Print ["***  MULTIPLE  COMMON  RATIOS?  ***"] 

]; 

ratio=Max[ common] ; 
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h2=l /ratio ; 

gl=(b01  cll  +  b02  cl2  h2''2)  /  (Sqrt  [lam]  *  (bOl  +  b02  h2)  )  ; 
g2=(b01  c21  +  b02  c22  h2''2)  /  (Sqrt  [lam]  *  (bOl  +  b02  h2)); 
{{gl,g2},h2} 

]  ; 

process [ v_, str_] : = 

Block [ { scale, ints , list , min, max, ten^, i , PI ] , 

Print ["Mean:  ",Mean[v]]; 

Print [ "Median :  " , Median [ v] ] ; 

Print ["Variance:  " , Variance [v] ] ; 

Print [ " Standard  Deviation :  " , StandardDeviation [ v] ] ; 

Print [ " Range :  " , SampleRange [ v] ] ; 

Print ["Minimum:  " ,min=Min[v] ] ; 

Print ["Maximum:  " ,max=Max[v] ] ; 
scale=0 ; 

■While  [Round [max*  10^ scale ]  —Round [min*  10'' scale] ■<10 ,  scale++  ]  ; 

While  [Rovind[max*10''scale]  -Rovind  [min*  10'' scale]  >99 ,  scale — ]  ; 

ints=Round [v*10''scale]  ; 

list={ }; 

min=Min [ints] ; 

max=M2uc[ints]  ; 

For [ i=min , i<=max , i++ , 

temp=N[{i/10''scale,Co^Jnt[ints,i] }  ]  ; 

If  [teit5>[  [2]  ]  !=0,AppendTo[list,  ten?>]  ] 

]  ; 

Pl=ListPlot [list, Plot Joined- >True, PlotLabel->str, 

PlotRange->[0, Automatic] ,Ticks->{ {N [min/ 10^ scale] , 

N[  (min+max)/(2*10''scale)  ]  ,N[max/10''scale] } ,  Automatic) , 
AxesOrigin->{Min [v] ,0]] 

]  ; 


NA=A;  NC=C; 

{b01,b02}=NB;  {r01,r02}=NR; 
xvec={x''4 ,  x‘'3 ,  x''2  ,x,  1 }  ; 

Blist=Rlist={ } ; 

{R,B}=base [MA,MC] ; 

ndist=NonnalDistribution [0 , stdev*NC [ [ind [ [1 ] ] / ind [[2]]  ]]  ]; 

For [i=l , i<=i ter , i++ , 

NC=C; 

NC [ [ind [ [ 1 ] ] , ind [ [2 ] ]  ] ] +=Random [ndist ] ; 

While[NC[[ind[[l]],'ind[[2]]  ]]<0  || 

NC[[ind[[l]] ,ind[[2]]  ] ]>2*C[ [ind[ [1] ] ,ind[ [2] ]  ]], 

Print["***  ERROR  -  OUT  OF  BOUNDS  ***"]; 

Print [Abs [NC[ [ind[ [1] ] ,ind[ [2] ]  ] ] 
-C[[ind[[l]],ind[[2]]  ]]  ]  / 

{stdev*C[ [ind[ [1] ] ,ind[ [2] ]  ]]),"  stand,  dev."]; 
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NC[ [ind[ [1] ] ,ind[ [2] ]  ] ]=C[ [ind[ [1] ] ,ind[ [2] ]  ]]+ 

Random [ndist] ; 

]  ; 

results=base [NA,NC] ; 

If [i==l, 

Rlist={ {results[ [1,1] ] } , {results[ [1,2] ] } } ; 

Blist={ results [ [2] ] }  , 

AppendTo [Rlis t [[!]], results [[1,1]]  ]  ; 

AppendTo [Rlist [ [2 ] ] , results [[1,2]]  ] ; 

^^pendTo [Blist, results [ [2] ]  ] 

]  ; 

]  ; 

Print["H[2]:  ",B]; 

process [Blist, "H[2]  Distribution"] ; 

Print [ "  " ] ; 

Print["G[l]:  ",R[[1]]  ]; 

process [Rlist [ [1] ] , "G[l]  Distribution"] ; 

Print [ "  " ] ; 

Print["G[2]:  ",R[[2]]  ]; 

process [Rlist [ [2] ] ,"G[2]  Distribution"] ; 

] 
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Statistics  from  Sensitivity  Analysis 
Fire  Allocation  &  Attrition  Method 
June  3, 1994 

The  following  data  was  obtained  by  the  sensitivity  analysis  routines  for  the  Fire 
Allocation  Method.  The  initial  data  for  these  routines  were; 

A={{.  1  ,.05}.{.075,.  1 }}  NB={50.25} 

C={{.075,  .05},{.025, .  1 }}  NR={1 00.40} 

Statistics  from  perturbing  the  C[1,1l  entry  by  a  random  draw  from  a  normal 
distribution  with  5%  standard  deviation. 

H[2]:  0.849133 
Mean:  0.850096 
Median:  0.849214 
Variance:  0.0000909999 
Standard  Deviation:  0.00953939 
Range:  0.0576209 
Minimum:  0.825638 
Maximum:  0.883259 

G[l] :  0.871234 
Mean:  0.870513 
Median:  0.871098 
Variance:  0.000288957 
Standard  Deviation:  0.0169987 
Range :  0 . 101812 
Minimum:  0.819872 
Maximum:  0.921684 

G[2]:  0.571781 
Mean:  0.573207 
Median:  0.571927 
Variance:  0.000336314 
Standard  Deviation:  0.0183389 
Range:  0.110249 
Minimum:  0.52304 
Maximum:  0.633289 

Statistics  from  simultaneously  perturbing  each  entry  in  the  first  row  of  C  by  a 
random  draw  from  a  normal  distribution  with  5%  standard  deviation. 

H[2]:  0.849138 
Mean:  0.849598 
Median:  0.848658 
Variance:  0.000134497 
Standard  Deviation:  0.0115973 
Range:  0.0714355 
Minimum:  0.82365 
Maximum:  0 . 895086 

\ 

G[l] :  0.871234 
Mean:  0.871896 
Median:  0.872109 
Variance:  0.000427619 
Standard  Deviation:  0.0206789 
Range:  0.121111 
Minimum:  0.305736 
Maximum:  0.926846 
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G  [2 ]  :  0.571731 
Mean:  0.572038 
Median:  0.570843 
y'ariance:  0.000496493 
Standard  Deviation:  0.0222321 
Range :  0.134293 
Minimum :  0.518575 
Maxiiaum:  0.6528  63 

Statistics  from  independently  perturbing  each  entry  in  the  first  row  of  C  by  a 
random  draw  from  a  normal  distribution  with  5%  standard  deviation. 

H[2]:  0.849133 
Mean:  0.849416 
Median:  0.848763 
Variance:  0.0000662001 
Standard  Deviation:  0.0092344 
Range:  0.0625471 
Minimum :  0.325587 
Maximum:  0.888134 

G[l] :  0.871234 
Mean:  0.871701 
Median:  0.871917 
Variance:  0.000278723 
Standard  Deviation:  0.0166951 
Range:  0.107946 
Minimum :  0.313867 
Maximum:  0.921313 


G[2i  :  0,571781 
Mean:  0.571909 
Median:  0.571052 
Variance:  0.000320908 
Standard  Deviation:  0.0179139 
Range:  0.118523 
Minimum:  0.522927 
Maximum :  0.64145 

Statistics  from  simultaneously  perturbing  each  entry  In  the  first  column  of  C 
by  a  random  draw  from  a  normal  distribution  with  5%  standard  deviation. 

H[2]:  0.849138 
Mean :  0.849248 
Median:  0.849159 
Variance:  0.0000196365 
Standard  Deviation:  0.00443131 
Range:  0.0264521 
Minimum:  0.339036 
Maximum:  0.865538 

G[l] :  0.871234 
Mean:  0.371637 
Median:  0.871154 
Variance:  0.000286002 
Standard  Deviation:  0.0169116 
Range:  0.0973707 
Minimum :  0.817552 
Maximum:  0.914923 
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G[2]:  0.571731 
Mean :  0.572005 
Median:  0 . 571793 

-6 

Vaiiance:  7.3137  10 

Standard  Deviation:  0.00270435 

Range:  0.0166747 

Mi  n  iruum :  0.566729 

Max  iiuuiu :  0.533403 

Statistics  from  independently  perturbing  each  entry  in  the  first  column  of  C 
by  a  random  draw  from  a  normal  distribution  with  5%  standard  deviation. 

H[2i:  0.849138 
Mean:  0.850224 
Median:  0.349267 
Variance:  0.000115734 
Standard  Deviation:  0.0107603 
Range:  0.0654364 
Minimum:  0.818917 
Ma-ximum :  0.884353 

G[l] :  0.871234 
Mean:  0.370493 
Median:  0.870342 
Variance:  0.000298903 
Standard  Deviation:  0.0172888 
Range:  0.100964 
Minimum :  0 . 3 1 6  9  7  3 
Maximum:  0.917942 

G[2] :  0.571731 

Mean :  0.573444 

Median:  0.572479 

Variance:  0.000563436 

Standard  Deviation:  0.0237368 

Range:  0.146793 

Minimum :  0.495818 

Maximum :  0.642611 
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Graphs  from  Sensitivity  Analysis 
Fire  Allocation  &  Attrition  Method 
June  3, 1994 


The  following  graphs  were  obtained  by  the  sensitivity  analysis  routines  for  the  Fire 
Allocation  Method,  The  initial  data  for  these  routines  were; 

A={{.  1 .  .05}, {.075, .  1 }}  NB={50,25} 

C={{.075,  .05},{.025, .  1 }}  NR={1 00,40} 

Graphs  from  perturbing  the  C[1,1]  entry  by  a  random  draw  from  a  normal 
distribution  with  5%  standard  deviation. 

H[2]  Distribution 


0.575 


0.63 


Graphs  from  simultaneously  perturbing  each  entry  in  the  first  row  of  C  by  a 
random  draw  from  a  normal  distribution  with  5%  standard  deviation. 

H  [  2  ]  Di s  tr ibutiori 


G[2]  Distribution 
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Graphs  from  independently  perturbing  each  entry  in  the  first  row  of  C  by  a 
random  draw  from  a  normal  distribution  with  5%  standard  deviation. 


H [2]  Distribution 
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Graphs  from  simultaneously  perturbing  each  entry  in  the  firet  column  of  C  by 
a  random  draw  from  a  normal  distribution  with  5%  standard  deviation. 

H[2]  Distribution 
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Graphs  from  independently  perturbing  each  entry  In  the  first  column  of  C  by  a 
random  draw  from  a  normal  distribution  with  5%  standard  deviation. 

H[2]  Distribution 
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This  is  just  a  sample  of  the  numerical  experiments  that  we  ran.  The  results  led 
to  no  substantial  conclusions,  except  perhaps  that  the  ATCAL  method  is  somewhat  less 
sensitive  to  changes  in  the  initial  data.  All  the  methods  produced  reasonable  results 
(although  not  consistent  results  across  the  methods).  Efforts  should  be  made  to  do  large- 
scale  comparisons  perhaps  using  a  set  of  COSAGE  replications  as  a  common  data  pool  for 
all  three  methods. 
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§5.  Conclusions. 

We  have  examined  three  methods  for  computing  weapon  scores  (values,  importances). 
Of  the  three,  the  fire  allocation  methods  are  the  most  general,  but  computationally  they 
are  the  most  difficult.  The  other  two  methods,  the  eigenvalue  or  potential/antipotential 
method  and  the  method  employed  in  ATCAL,  are  both  in  use  in  existing  combat  simula¬ 
tions. 

Among  the  results  we  have  obtained  is  the  proof  of  the  existence  and  uniqueness 
for  solutions  to  the  ATCAL  importance  equations,  under  some  mild  conditions.  This 
leads  to  a  robust  iterative  scheme  for  solving  these  importance  equations.  While  the 
ATCAL  equations  are  somewhat  ad  hoc  in  nature,  they  do  give  reasonable  scores  for 
various  weapons  that  are  generally  consistent  with  the  other  methods.  (Our  numerical 
experiments  were  limited  to  a  small  number  of  weapon  types  on  each  side  (two  or  three) 
which  limits  our  conclusions  involving  comparisons  across  methods.) 

Our  most  striking  observation  is  that  all  three  methods  seem  to  be  based  on  similar 
mathematical  phenomenon.  The  potential/ antipotential  method  and  the  fire  allocation 
methods  both  result  in  a  kind  of  “simultaneous  eigenvalue  problem”  -  a  real  simultaneous 
eigenvalue  problem  in  the  first  case  and  a  non-linear  analogue  in  the  second  case.  The 
potential/antipotential  method  and  the  method  in  ATCAL  rely  on  virtually  identical  proofs 
for  existence  and  uniqueness,  both  are  centered  around  a  common  fixed  point  theorem. 
Thus  ATCAL  yields  a  different  non-linear  generalization  of  the  linear  eigenvalue  method. 

The  good  news  is  that  these  similarities  give  us  some  assurance,  albeit  weak  assurance, 
that  all  three  methods  will  yield  generally  consistent  results.  This  is  confirmed  by  our 
relatively  limited  numerical  experiments.  Nevertheless  the  methods  do  not  give  identical 
results  and  caution  is  always  in  order  whenever  comparisons  are  being  made,  because 
results  obtained  from  one  method  are  not  directly  comparable  to  results  obtained  from 
another. 

The  bad  news  is  that  we  can  also  expect  the  methods  to  suffer  from  similar  qualitative 
defects  and  degeneracies.  One  will  probably  be  able  to  create  examples  where  all  the 
methods  fail  to  give  “reasonable”  answers,  if  you  can  create  an  example  for  any  one  of  the 
methods.  More  specifically,  many  of  the  criticisms  leveled  at  the  eigenvalue  approach,  in 
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all  likelihood,  apply  (with  appropriate  modification)  to  the  other  two  methods. 

The  author  would  again  like  to  express  his  thanks  to  the  U.S.  Army  Concepts  Analysis 
Agency  and  the  Army  Research  Office  for  their  support  of  this  research. 
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